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A MATHEMATICAL ANALYSIS OF OXYGEN 
RESPIRATION IN MAN 


ARTHUR B. CuILton,* DELBERT S. BARTH,{ 


AND RALPH W. STACY 


DEPARTMENTS OF Puysics, PHysics, AND PHYSIOLOGY, RESPECTIVELY 
OxI0 STATE UNIVERSITY 


By means of physicomathematical models similar to those used by A. B. Chilton and 
R. W. Stacy (1952) formulas are obtained for the purpose of analyzing oxygen respiration in 
man under conditions of metabolic equilibrium. The formulas are applied to conditions in- 
volving a variety of physiological assumptions. Calculated quantities are shown to be in agree- 
ment with experimental values. A quantitative prediction is made of the fluctuation in alveolar 
oxygen tension over a single respiratory cycle. Finally, calculations are completed for allalveolar 
gases; and the over-all results are shown to be self-consistent to the degree expected in view of 
the simplifying approximations made. 


I. Introduction. This paper is a companion to a similar one on carbon 
dioxide respiration written by Chilton and Stacy (loc. cit.). For a more 
complete statement of the background for a mathematical approach to 
respiration, as well as a compilation of the assumptions required to give 
a model of the respiratory system capable of analysis, a study of the 
original paper is recommended. 

In fashion similar to that of the paper on carbon dioxide, we propose to 
analyze the kinetics of oxygen respiration in the lungs under conditions of 
metabolic equilibrium. The first step will be to solve a basic problem as- 
suming constant alveolar volume; the next will be to generalize to cases 
of inspiration and expiration. The cyclic changes in oxygen tension in the 
alveoli will be determined; and finally the effect on the results of changes 
in the various factors assumed originally. 

A supplementary section is included which contains certain minor 
points of discussion, for the ultimate purpose of providing an over-all 
check upon our theoretical results. 

* Commander, Civil Engineer Corps, U.S. Navy, attending Ohio State University under 
the sponsorship of the Office of Naval Research and the U.S. Naval Postgraduate School. 


+ Captain, Chemical Corps, U.S. Army, attending Ohio State University under the super- 
vision of the U.S. Naval Postgraduate School. 
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II. Analysis of oxygen flow in model (volume constant). The initial ap- 
proach to the problem is the same as in the carbon dioxide paper, so that 
a brief outline is sufficient up to the point that a change in approach is 
necessary. The symbols used will be the same, wherever possible. 

Figure 1 depicts the model giving in schematic fashion the physical re- 
lationships considered to exist in the lungs. Figure 2 gives the simplified 
model used for most of the analysis. 

Let p be the partial pressure of CO, in the alveolar space having vol- 
ume V. Let M be the molar quantity of oxygen in space V. The gas con- 


| 


Frcure 1. Schematic diagram of the physical relations in the lungs. Space U represents 
that portion of the lungs which is ventilated but not perfused with capillaries; V is that portion 
which is both ventilated and perfused; Y is perfused but not ventilated; Z is neither perfused 
nor ventilated. The process of breathing is depicted as a movement of a piston (hatched) within 
a cylinder. (Same as Fig. 1, Chilton and Stacy, Joc. cit.) 


tent is R, and T is the absolute temperature. Then from the ideal gas law, 
it may be stated that 


pV = MRT. (1) 


Further, the following symbols are defined: v is the velocity of blood 
flow through the capillary channel, L is the length of the channel, and B is 
the cross-sectional area of the channel. Subscripts ; and 2 refer to the con- 
ditions of the blood at the entrance and exit, respectively, of the capillary 
channel. Let c be the concentration of oxygen in the blood: é refers to the 
average value of c throughout the capillary channel at any instant of 


time. The subscript 9 refers to known initial conditions, when time (é) is 
zero. 
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From Figure 2, it is easy to find the rate of increase of oxygen in the 
alveolar volume: 
dM dé 


F =Bv(¢,— Co) —BL —. (2) 


(Note: This equation is not precisely correct, as another source of oxygen 
exists in the system. The nature of this factor and the degree of approxi- 
mation occasioned by ignoring it are more easily explained after the re- 
sults of the more approximate solution are obtained. Such discussion is 
therefore relegated to the final section of this paper.) 


FOR CALCULATIONS 
USING CONSTANT 


VOLUME, THIS 
PISTON IS ASSUMED 
STATIONARY. 


Ficure 2. Simplified diagram of physical relations. (Same as Fig. 2, Chilton and Stacy, 
loc. cit.) 


In order to solve this equation, we must first find dé/dt. Figure 3 gives 
the O.-dissociation curve for the blood, so that one can relate oxygen con- 
centration to oxygen tension in the blood as it passes through the capil- 
lary channel. From this figure we note that if the oxygen tension of blood 
is 80 mm. Hg or over (as it is normally for blood leaving the lungs), the 
concentration of oxygen in the blood is essentially constant. Thus we take 
c2 to be known and to be constant over the respiratory cycle. The assump- 
tion of metabolic equilibrium, of course, prescribes that c; is also constant. 

To determine é, we should know the variation of c with distance along 
the capillary channel. It is well known that even under conditions of severe 
exercise the blood leaves the lungs with a full complement of O, (Nims, 
1949, p. 803). This leads us to a realization that, under normal conditions, 
the variation of oxygen concentration in the alveolar capillary channel is 
approximately as indicated in Figure 4. 
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We see that, to a rough degree of approximation, 


€ =a" (3) 
and thus 

dc. des _ 4) 

IPP Cae ( 


It should be noted that equation (4) is even more correct than equation 
(3), for the effect of cyclic variations in oxygen tension in the alveoli is 
very small in a precise determination of é (see Fig. 4). For the sake of the 
reader who may be familiar with the concept of “effective volume” ex- 
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Ficure 3. Diagram of oxygen dissociation curves corresponding to pCOz levels found in 
arterial and venous blood, and the physiological Op» dissociation curve. 


plained by Chilton and Stacy (Joc. cit.) for carbon dioxide respiration 
analysis, we are now saying essentially that for the oxygen case the effec- 
tive volume of the lungs differs little from the actual volume. (We have 
computed that effective volume is only on the order of 1% above actual 
volume for oxygen, although about 12% more for carbon dioxide.) 

On the basis of equation (4), then, we may drop the last term in equa- 
tion (2), so that our differential equation becomes 


OM =Bo(c— a). (5) 
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If we replace M by its value from equation (1), we obtain 


where 
p = RT - Bo (co — ¢). Gi) 


Having a given initial condition for p and considering that V is constant, 
we find that 
p = Po _— 7 ¢ (8) 


III. Analysis of oxygen tension in model (variable volume). With the 
same qualifications that were made in the carbon dioxide case (though 
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Ficure 4. Variation of O2concentration of blood in alveolar capillaries at two different times 
in the respiratory cycle. The horizontal scale (X) represents distance along the capillary. 


they are still less important here), we generalize to the case of variable 
volume by assuming that equation (6) is still valid when V is dependent 
upon time. 

(A) Variable volume without dilution, as during expiration: Suppose 


ie — Vo —— bs (9) 
Then substitution in equation (6) gives 


Sa Renee aie (10) 
dt Vo—mt 


Having a given initial condition for p, we find that 


p= mt, in [1-7]. (11) 
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(B) Variable volume with dilution, as during inspiration: Suppose 
V=Votm. (12) 


Then, by exactly the same line of reasoning as was used for the carbon 
dioxide problem, we obtain the following differential equation 


dp~ = j Some? 2 
dt Votmi NES is 
m 
(13) 
(»'-4)-» 
pee 
m 


where p’ (considered a constant in this equation) is the pO: of the dilut- 
ing gas. 
Having a given initial value of p, we find then that 


= a oa) + i mie es he 


The value of p’ depends upon the source of the inspired gas. If the in- 
spired gas is saturated air, then p’ is the partial pressure of saturated- 
atmospheric oxygen (149 mm. Hg). During the initial stage of inspira- 
tion, however, the inspired gas has the average composition of the previ- 
ously expired gas which was retained in the dead space. If this gas repre- 
sents that expired from the alveolar volume during the period, t; to fs, 


P= fi mott in -F a, (15) 
which can be reduced to the following form: 
sults ms dm 
rebate 
RA seal grey by 


In equations (15) and (16), the subscript » refers to the conditions at 
the beginning of expiration, whereas in (14) the subscript o refers of course 
to the conditions at the beginning of an inspiration phase. 

IV. Application of derived formulas. In order to apply the formulas 
derived above, it is necessary to list the values of the physiological con- 
stants entering into the equations. Where applicable they correspond to 
the values assumed for the carbon dioxide problem: 
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R is 62,400 mm. Hg — cc./moles — °K; T is constant at 310° K; 
Bo is taken to be 83.3 cc./sec. under normal resting conditions, and 250 
cc./sec. during exercise. 

The oxygen content of the blood is taken as 13 vol. % as it enters the 
lungs (normal conditions) and 19.5 vol. % as it leaves. Under exercise 
conditions the oxygen content of the entering blood is quite low, and taken 
as 3.5 vol. %. This means that in units of moles/cc. the value of c; (nor- 
mal) is 5.8 X 10-, the value of ¢ (exercise) is 1.56 X 10-°, and the value 
OL 1s. 8.7 < 10-*. 

Then the value of u (normal) is 4,675, and of u (exercise) is 34,500. 
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Ficure 5. Cyclic variation of alveolar pO2 and alveolar volume in “‘normal” respiration 


Consider lung volume at inspiration as 3440 cc. and at expiration 2890 
cc., giving a tidal volume of 550 cc. If we take dead space volume as 165 
cc., then V varies between 2725 cc. and 3275 cc. In exercise we assume V 
varies from 3500 cc. to 2000 cc., with a dead space of 250 cc. 

The length of a single respiratory cycle is taken to be 3.5 seconds, but 
in the case of exercise a 2-second cycle is assumed. 

Figure 5 shows the respiratory pattern consistent with the data. For 
computational purposes, the breathing pattern is approximated rather 
closely by a series of straight lines. 

In the phase a — 8, the expiration phase, 


V = 3275 — 3671, bdsig 
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and by use of equation (11), we can find 


Pepa aoe (18) 


In phase b — c (the respiratory phase) V is 2725, and by equation (8) 
we find 
pe = fr — 0.86. (19) 


Phase c — d, the first part of inspiration, is that part wherein dead 
space air, accumulated during the latter part of expiration, is re-inspired. 
The time of this phase is 0.45 seconds. The value of p’, the pO, of this 
re-inspired gas, is given by equation (16) and is found to be 

b' = ~,—1.93. (20) 

In phase c — d, 

V = 3275 + 367 4, (21) 
and by use of equation (14), we can find 


ba = .057 pat .943 p. — 0.84. (22) 


In phase d — a, the latter part of inspiration, 
V = 2890 + 367 ¢t. (23) 
The value of p’ is determined for atmospheric air, fully saturated with 
water vapor at body temperature, since saturation is almost instantane- 


ous for air as it passes through the trachea to the alveoli. By use of equa- 
tion (14), p’ being taken as 149 mm. Hg, we obtain 


pa = .882 pa + 16.03. (24) 


Equations (18), (19), (22), and (24) constitute a set of four equations 
in four unknowns which can be solved simultaneously to give 


ba = 107.3 

pe = 104.95 

pe = 104.1 Sod 
ba = 103.4. 


These values are plotted in Figure 5, below the respiration pattern curve. 
Note that the compound curve from , to a is practically a straight line; 
this indicates that under normal circumstances, a single formula may give 
the results for pO: from a — d. In fact, as the reader may verify, the use 
of equation (9), with the mean alveolar volume in phases, a — d, will give 
almost the same predicted values as (25). 
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Since we now have formulas for instantaneous value of pOnz, it is possible 
by an integration process to determine the flow of oxygen in and out of the 
alveolar volume in the breathing cycle, and thus compute the net rate of 
oxygen consumption of the body. This approach is desirable only as a 
check on the accuracy of our previous computation, for a simpler method 
of obtaining an accurate value to the oxygen consumption rate is available 
[see eq. (5)]. 

We have assumed that the normal individual’s blood enters the lungs 
with a concentration of oxygen of 13 vol. % and leaves with 193 vol. %. 
This represents a difference of 63 vol. %. If 5000 cc. of blood pass through 


TABLE I 


THEORETICAL ANALYSIS OF THE INFLUENCE OF CHANGES IN SOME PHYSIOLOGICAL 
FACTORS IN O2 CONSUMPTION AND ALVEOLAR pO, 


CHANGE IN PHysrIoLocicaL Factors CHANGE In CoMPUTED ReEsuLTS 
anes Os Median Cyclic 
Altered Factor Amount Consumption Alv. p02 Variation 
Rate (mm. Hg) | (mm. Hg) 
1. Arterial O2 Vol. % + 0.5 + 7.7% — 3.4 + 0.3 
2 Venous O2 Vol. % — 0.5 + 7.7% — 3.4 + 0.3 
se Total Blood Flow Rate + 7.7% + 7.7% — 3.4 + 0.3 
4. Av. Alveolar Volume — 8.3% 0 — 0.3 + 0.4 
ss Breathing Rate +75 %G 0 +19.2 — 1.9 
6. Tidal Volume —-9 G 0 — 6.6 + 0.1 
he Dead Space —-9 GWG 0 + 1.7 0 
8. Moderate Exercise 2529 (it |-846... ek): GAD) (Gti): SRE a ee +11.2 
9. Inact. Part of Alv. Vol. 9G 0 — 2.4*| + 0.2* 
10. Enrich pO: in Air Breathed +6% 0 + 9.2 0 


* Space V only. See text. 


the alveoli per minute, the amount of oxygen emitted is (5000/100) X 65, 
or 325 cc./min. (STP). This figure can be verified by the more tedious 
method outlined in the preceding paragraph. 

The median pO, is easily seen to be 104.85 mm. Hg, and the cyclic 
variation is 3.9 mm. Hg. 

The previous calculations and results have been based on what we call 
the “normal” individual’s characteristics while at rest. Since our assump- 
tions are to a certain extent arbitrary, it is quite desirable to determine 
the differences to be expected in the computed results if one characteristic 
is changed for each separate calculation. Table I lists the variations in 
data used and results obtained. In the table, data are given relative to the 
basic solution for “normal” conditions. With the exception of Case 9, the 
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calculations follow almost precisely the pattern presented above. (Note: 
Care should be taken not to use Table I when the pO; is less than about 
80 mm. Hg, since the assumption of constant ¢2 is no longer valid under 
such circumstances.) 

This case assumes that 9 per cent of the alveolar volume is considered 
as inactive in the calculations. Thus our model here is more complex than 
that given in Figure 2, since it includes both spaces U and V of Figure 1. 
The special features to be considered in the calculations for this case are 
noted in the carbon dioxide paper. Results obtained for this case are 


pa = 104.4 ) 

ps = 101.8 

pe = 100.9 (26) 
pa = 100.3 


fa = fo = fe = 135.6 


(f refers to pO» for space U, as p does for space V). 

V. Discussion of results. (A) The oxygen consumption rate for the nor- 
mal case, as indicated above, is 325 cc./min. This is roughly within the 
normal limits prescribed by physiological literature. Taken in conjunction 
with the CO, output rate of 272 cc./min. (Chilton and Stacy, loc. cit.), 
this figure represents a value of RQ (Respiratory Quotient) of 0.84, which 
is about average for human beings. Such an exact prediction of theory is 
somewhat fortuitous, but any reasonable set of basic physiological as- 
sumptions will provide an approximately valid figure for RQ. 

(B) The prediction of a cyclic variation in pO, in the active alveolar 
space (amounting to a little less than 4 mm. Hg), just as for the previous 
prediction of pCO, variation, has not yet been adequately measured ex- 
perimentally. Such a variation becomes much less for forced rapid breath- 
ing (Case 5), but much greater under exercise conditions (Case 8). 

(C) It may be noted that the rate of oxygen consumption is essentially 
constant over the cycle, since none of the factors entering into the com- 
putation of oxygen consumption rate changes within the period of a single 
respiratory cycle. 

Furthermore, since oxygen consumption is dependent on very few 
variables [equation (5)], many variations in the physiological constants 
have no effect on oxygen consumption rate. This is in some contrast to the 
carbon dioxide behavior. For example an increase in breathing rate (forced 
breathing) has a strong effect on the rate of carbon dioxide elimination, 
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but practically none on oxygen uptake. (As usual, our computations do 
not take into consideration any secondary effects via reflex pathways or 
chemical action on the respiratory center.) 

(D) The theoretical prediction of 104.85 mm. Hg for average alveolar 
pOz is fairly close to the traditional experimental value of 100 mm. Hg. It 
must be stated moreover that recent experiments at this institution (Ohio 
State University) with use of the mass spectrograph (Stacy e¢ al., 1948) 
indicate that the traditionally accepted value for alveolar pO, is indeed 
a little low. This is considered as excellent evidence of the validity of our 
approach to the oxygen respiration problem. 

(E) A knowledge of the cyclic variations in the gas concentrations of 
the blood as it enters and leaves the lungs enables us to determine a cyclic 
variation in instantaneous RQ, defined as the ratio of instantaneous rate 
of carbon dioxide removal from the blood to the rate of oxygen uptake 
by the blood. If we consider that carbon dioxide output rate is proportion- 
al to the difference in concentration (or volume per cent) of the incoming 
and outgoing blood of the lungs, the instantaneous rates are seen to be 
(see Chilton and Stacy, 1952, p. 14 and Fig. 5): 


At instant a, CO: output rate = 325 cc./min. 


At instant 6, CO, output rate = 259 cc./min. (27) 
At instant c, CO: output rate = 238 cc./min. 


At instant d, CO2 output rate = 223 cc./min. ) 


Since the O, consumption rate is constant at 325 cc./min., as we have 
seen, the instantaneous values of RQ become: 


At instant a, inst RQ = 332 = 1.000 
At instant 6, inst RQ = 332 = .797 (28) 
23 


At instant c, inst RQ = 238 = .733. 


It must not be considered that these instantaneous values of RQ have 
the general metabolic significance that the average RQ has. 

VI. Supplementary comments. Throughout the foregoing analysis a fac- 
tor has been ignored, and such an action requires justification, since its 
importance, or lack of same, is not immediately evident. This factor is 
explained as follows: 

Consider the system with constant alveolar volume. Equation (2) is 
the mathematical expression for conservation of mass, with sources (or 
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sinks) of oxygen at points 1 and 2 in the alveolar capillary channel fF ig. 
2). We have neglected the possibility that an oxygen source may exist at 
the opening to the dead space. (We do not refer to diffusion between vol- 
ume V and the dead space—this is considered inconsequential.) We note 
that with an RQ of about 0.84, the system is losing a greater volume of 
oxygen to the blood than the volume of carbon dioxide it receives from 
the blood. If pressure in volume V remains at the atmospheric value, a 
certain amount of dead space gas passes into the alveolar space to make 
up for the net loss of gas in the lungs. 

Since the instantaneous RQ in the respiratory cycle.is practically al- 
ways different from unity, this effect should in theory be considered in all 
phases of the cycle. It becomes evident, however, that instantaneous RQ 
varies over the period of the cycle, so that precise equations become rather 
elaborate and exact calculations become quite tedious. Such calculations 
have been carried out by the authors, but the results obtained differed 
little from those given in the body of this paper. The corrections may also 
be provided rather accurately by the following approximate procedure, in- 
volving the use of Table I. 

Consider the constant volume case. The rate of loss of oxygen to the 
blood is 325 cc./min. (STP); and this is equivalent to the average value 
of dM /dt throughout the respiratory cycle, according to our basic analysis 
above. The rate of carbon dioxide accumulation from the blood is 272 
cc./min. (STP), and the resulting difference is 53 cc./min. (STP). Thus 
gas from the dead space flows into the active alveoli at about this rate, 
or at least during inspiration and expiration the rate of gas inspiration or 
expiration is modified by this amount. This amount of gas is readily seen 
to have a pO: of about 105 mm. Hg during phases a — d, and 149 mm. 
Hg during the phase d — a. This gives an average pO». of about 118 
mm. Hg over the whole system. 

Thus additional oxygen enters the alveolar volume at the rate of 
(118/760) X 53 cc./min., or 8.23 cc./min. This corresponds to an average 
increased rate of 8.23/325, or 2.53%. Roughly, then, equation (5) may 
be corrected by increasing the term dM/di by a factor of 1.0253. Such a 
factor can then be transferred to the other side of the equation, and the 
given development can be followed, with the use of py’ replacing un, where 


qetiog 
7.0253. (29) 


This simple substitution can be carried through all conditions specified 
of inspiration and expiration, so that the corrected calculation is precisely 
like the original calculation, with y decreased by about 23%. 
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The results from such a calculation can be obtained from Table I. The 
first three cases of the table involve altered factors which essentially in- 
volve nothing more in the calculation than an increase of 7.7% in p. We 
may presume that our corrections should be 2.5/7.7 of those listed, but 
with the opposite sign. 

This means that the median O; level in the alveoli is higher than that 
originally computed by about 1.1 mm. Hg, and the cyclic variation in 
pOz is decreased by about 0.1 mm. Hg. No correction in the rate of oxygen 
consumption is valid, since the rate must be based on the original u rather 
than the corrected py’. 

The corrections are thus seen to be too small to make a more elaborate 
analysis in this paper desirable. 

It should be noted that in the case of carbon dioxide, such a correction 
is still less necessary. The average pCOy of dead space gas is about 25 mm. 
Hg, so that the average rate of influx of CO, for the dead space is about 
(25/760) X 53, or 1.76 cc./min. This amounts to only 0.65% of the aver- 
age rate of influx of CO, flow from the blood to the alveoli. 

Finally, we may roughly compute the pNp of the alveolar air. Since the 
influx of gas into the lungs to replace the net loss of 53 cc./min. (STP) is 
largely nitrogen, it can be seen that the alveolar gas has an appreciably 
higher concentration of nitrogen than does the saturated atmospheric air. 

The nitrogen lost in exhalation per cycle is accomplished during the 
first 1.05 seconds of phase a — 0; and is, in cubic centimeters, 1.05 X 
Pav/760 X 367, where Pav is the average alveoli pNo. That nitrogen emit- 
ted from the alveolar space in the latter 0.45 seconds of phase a — 0 is 
retained in the dead space and re-inhaled. It may be disregarded. The 
nitrogen gained by inhalation, during the latter 1.05 seconds of phase 
d —a, is 1.05 X p’/760 X 367, where p’ is the pNe of saturated at- 
mospheric air (= 564 mm. Hg). The net loss in a cycle is 1.05 X 
(Pav — p’) X 3677/60 cubic centimeters. This is balanced by the nitro- 
gen gain in the gas influx to make up for the 53 cc. (STP) per minute, or 
3.5 cc. (body temperature) per cycle lost in the blood gas exchange. This 
amounts to 3.5 X ’/760 cc. per cycle. 

Equating these two values gives 

385 3.5 X564 


oe ee TE 


30 
760 eee 


This gives us the value of Pav as 569.1 mm. Hg. 
The accuracy of computation on all gases may be checked, to the ex- 
tent permitted by our approximations, by adding the values of gas pres- 
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sures computed at each instant of the alveolar cycle. Table IT gives sum- 
mations at instants a, b, c, and d in the cycle. The results are sufficiently 
close to 760 mm. Hg to give some confidence in the self-consistency of 
our calculations. Since we have neglected small cyclic variations to be ex- 
pected in pNb, the variations in the totals are not considered significant 
of appreciable error in pO: or pCO» cyclic variations predicted. 

Table II also is a summary of our theoretical results for all alveolar gas 
tensions in the “normal” individual. 


TABLE II 


VARIATIONS IN ALVEOLAR CONDITIONS DURING 
A SINGLE RESPIRATORY CYCLE 


INSTANT IN CYCLE 


Gas Tension a b c d 
D CO a saenerer ae Para 34.9 36.7 my ee: 37.8 
pOz (basic comput.)...... 107.3 105.0 104.1 103.4 
Correction to pO2...... + 1.1 + 1.1 + 1.1 + 1.1 
DIN xt here eas oe PS 569.1 569.1 569.1 569.1 
PHeO cern d seat eee 47.0 47.0 47.0 47.0 
‘LOtaley aera: shorten: 759.4 758.9 758.6 758.4 
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Two first-order, non-linear differential equations are presented to describe the interaction 
of the thyroid and pituitary glands and to explain some of the phenomena of periodic relapsing 
catatonia. Topographical study of these equations shows a mechanism in which either random 
fluctuations of hormone levels or system stability may occur. The variations in hormone con- 
centrations are assumed to be caused by system malfunction, while stability is shown to result 
from the administration of exogenous thyroid which, if supplied at an optimum rate, suppresses 
the output of both glands. The system, under treatment, shows a zero level of thyrotropin and 
a thyroid concentration which is larger than the equilibrium level in the untreated system. The 
theory is consistent with what is known of the course and proper treatment of periodic relapsing 
catatonia. 


Periodic relapsing catatonia, first described and successfully treated by 
R. Gjessing (1939), is a type of mental disorder marked by fairly regular 
variations in the patient’s behavior. These variations persist, usually until 
death, unless the patient is given appropriate amounts of thyroid extract 
or thyroxin; and they return if the daily intake of thyroid falls below a 
level which is critical for the patient. Shock therapy may relieve the 
symptoms temporarily, but often is without effect and is never of lasting 
benefit. 

Gjessing (loc. cit.), Danziger et al. (1948), A. G. Gornall et al. (1953), 
and G. Mall (1952) observed that the clinical condition is correlated with 
fluctuations in the basal metabolic rate (BMR). In very elegant experi- 
ments, Gjessing found regular changes in the nitrogen balance, which we 
believe reflect changes in the BMR. 

In explaining these findings, we postulate the following concepts of 
thyroid-pituitary interactions: The pituitary makes a hormone, thyrot- 
ropin, which stimulates the thyroid gland to make thyroid hormone, 
which in turn inhibits the production of thyrotropin. If the pituitary is 
removed, the production of thyroid hormone falls to zero, while if the 
thyroid gland is removed, the production of thyrotropin proceeds un- 


* Read before the American Chemical Society, Chicago, Tll., September 11, 1953. 
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checked at a rate which we assume is constant. If a sufficient quantity of 
thyroid extract is administered by a physician, the production of thyrot- 
ropin is driven to zero. The two hormones are assumed to stimulate and 
inhibit in accordance with the Langmuir adsorption isotherm (Langmuir, 
1916) and are assumed to be excreted or otherwise lost at rates which are 
proportional to the hormone concentrations. These ideas may be expressed 
mathematically by the following non-linear differential equations: 


db kimt 
dt itm 


dr _ — kon@ vhs 
ie 7 


— b0=P(8, x) () 
gr =Q (8, wr) (2) 


PSU 20.5 


where 


6 is the level of thyroid hormone in the system at any time /, 

m is the level of thyrotropin-in the system at any time /, 

c is the rate of production of thyrotropin in the absence of thyroid inhibition, 

k; is a constant equal to the theoretical maximum production rate of the thyroid 
gland, 

kz is a constant assumed to be’greater than c so that the production of thyrotro- 
pin may be zero for sufficiently large @, 

m and n are the constants in the Langmuir adsorption equations, and 

b and g are the loss constants. 


These equations state that the rate of change of either hormone concen- 
tration is equal to the difference in the rate of production and the rate of 
loss of the hormone. As the production rates of either hormone depend 
upon the level of the other hormone, (1) and (2) describe a feed back sys- 
tem which will be shown to be a stable regulating mechanism whose ap- 
parent function is to maintain the equilibrium level of the two hormones. 

The action of the regulating system governed by equations (1) and (2) 
may be shown conveniently by plotting the locus of successive states in 
the 6 — 7 plane. Eliminating time from (1) and (2) yields 


dx _Q(8, x) ' 
dé. Pie,«)* (3) 


which describes the family of integral curves or the trajectories of the rep- 
resentative point (9, ) in the 6 — z plane. A topographical solution of 
(3) for any initial point (0, xo) may be obtained by the method of iso- 
clines (Minorsky, 1947). A simplified analytical procedure employing the 
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concept of isoclines will be employed to show the general shape of any 
trajectory. The two limiting isoclines, dx/d@ = 0 and dr/d@ = ~, may 


be obtained by setting (3) equal to zero and infinity respectively, yielding 
for dx/d@ = 0 the curve, 


QO(6, 7) = 0, (4) 
and for dr/d@ = = the curve, 
P(é, r) =0. (S) 


As shown in Figure 1, the limiting isoclinal curves divide the 6 — 7 plane 


0 8 —¢/n(ka-c) 


Ficure 1. The 6-x plane diagram and typical trajectory for equation (3) 


into four regions, which by analysis of (1) and (2) may be characterized 
by the following: 


Region A oe Le cs 0; > 0 (6) 
Region B a>; ts sea (7) 
Region C <0; aE <0} 0 (8) 
Region D a <0; ane ey, peck (9) 


The shape of any trajectory defined by (3) and the initial point (0, To) is 
now clearly defined by the required values of the trajectory slopes on the 
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limiting isoclines and by the signs of these slopes in the four intermediate 
regions. Clearly, all trajectories must be clockwise spirals terminating at the 
singular point (6’, 7’) which is the equilibrium condition for the system. 
In the time domain, the spiral trajectories would appear as damped oscil- 
lations of 6 and about the equilibrium levels 6’ and x’. Figure 1 includes 
a typical trajectory as well as the trajectory directions on the limiting iso- 
clines and in the four intermediate regions of the @ — 7z plane. 

The action of the thyroid-pituitary regulating system as illustrated by 
this analysis can be summarized as follows: If, due to any disturbance, the 
hormone concentrations are driven away from the equilibrium point, the 
regulating action will return the system to equilibrium by means of rapidly 
damped oscillations of both hormone levels. 

In periodic relapsing catatonia, observed variations in the level of 
thyroid hormone are sustained with irregular amplitudes. This suggests a 
recurring malfunction of some part of the system which causes the hor- 
mone levels to be driven away from equilibrium. The response of the sys- 
tem of equations (1) and (2) to random disturbances would, if the equa- 
tions were valid in the intervals between disturbances, consist of trains 
of damped oscillations of the hormone levels. Although the source of in- 
stability is not known, stabilization of the thyroid level can be achieved 
by rendering the regulating system partially or wholly inoperative by 
substituting externally added thyroid for all or most of the internally pro- 
duced hormone. If thyroid is added to the system by daily doses of thy- 
roxin or thyroid extract, the equilibrium level of thyrotropin will be re- 
duced or suppressed entirely while the equilibrium level of thyroid hor- 
mone will be increased. These general effects have been partially confirmed 
by clinical observations and will be shown to result from mathematical 
analysis of the thyroid-pituitary system modified by an input of exogenous 
thyroid. 

Assume the system of equations (1) and (2) is modified by the daily ad- 
dition of equal amounts of thyroid extract. As the system time constants 
are long compared to the 24-hour period between doses, the intake of exog- 
enous thyroid extract can be assumed to be equivalent to a second source 
of thyroid whose rate of production is constant at a value taken, for con- 
venience, at DK where K is a constant and 3 is the loss constant in (1). 
As the rate of change of the thyroid level depends upon the equivalent 


rate of production, equation (1) will be modified by the factor 6K and can 
be written as 


dé kumar 


Wi Teas DIE ore: (10) 
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Equation (2) does not change, but equation (3) will contain the factor bK 
and becomes 


dr Q(8,7) 


d0 P(@—K,z)~ ae 


The trajectories of the modified system are now defined by (11) and can 
be sketched in a manner similar to that previously employed. In the 
6 — mr plane diagram, the isocline dx/d#@ = 0 is unchanged as Q(0, 77) does 
not contain the parameter K, whereas the isocline dr/d@ = = is shifted 
to the right by an amount K as P(@ — K, x) in (11) differs from P(6, 7) 
in (3) only in that 6 — K is substituted for 0. It follows, therefore, that 


TT 


Ficure 2. The 6-7 plane diagram and typical trajectory for equation (11) with K 
= K’ =c/n(ke — ¢). 


the equilibrium point (6’, x’) for the modified system is shifted so that the 
equilibrium level of thyrotropin is reduced while the equilibrium thyroid 
level is increased. 

Consider the critical case where K = K’ = c/n(ke — c), the value of the 
9-intercept of the unchanged isocline dr/d0 = 0. The @ — 7 plane diagram 
for this case is sketched in Figure 2. Here the singular point (K’, 0) lies - 
on the 6-axis and at equilibrium, the output of both glands is zero while 
the system thyroid level is K’ and is supplied entirely from the external 
source. At equilibrium, therefore, the daily intake of thyroid just balances 
the loss and neither gland is required to deliver hormones to the system. 
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The response of this modified system with K = K’ is now limited to 
possible momentary output of either or both glands or a change in the 
loss constant 6. A disturbance of the first type would produce an initial 
point in the part of region B where 0 > K ‘ or in region C, and hence would 
be followed by a limited spiral trajectory continuing until the representa- 
tive point reaches the 6-axis. At this point, (01, 0), equations (2) and (10) 
would degenerate to 

a=0 (12) 
and 
oo + b6= dK’. (13) 


The representative point would, therefore, move along the 6-axis toward 
the equilibrium point (K’, 0) in accordance with the simple exponential 
time function which is the solution of the linear equation (13). A change 
in the loss constant 6 would modify the slope and position of the isocline 
dr/d@ = © and would cause a shift in the equilibrium point. The trajec- 
tory produced by an increase in 6 would be an exponential time decrease 
of @ toward the new equilibrium point with 7 remaining zero; and that 
produced by a decrease in b would be a small spiral originating in region 
A and proceeding to the new equilibrium point on the @-axis. A change in 
the loss constant g would produce no response as the equilibrium level of 
am is zero. 

As clinical observations show that the modification of the system as de- 
scribed here stabilizes the thyroid level without significant subsequent 
fluctuations, it is probable that the actual disturbing mechanism found in 
periodic relapsing catatonia is of the type which would produce no re- 
sponse in the modified system; namely, a decrease in the loss constant g 
or a momentary decrease in the output of either gland. 

The transient effects following the initiation of treatment or a change 
in the magnitude of the daily dose of thyroid extract remain to be dis- 
cussed, A treatment procedure in which a large initial dose or a series of 
large doses is administered, followed by smaller daily amounts to maintain 
the level established, would result in a large overshoot of the total thyroid 
level which could easily reach dangerous magnitudes. On the other hand, 
smaller doses, gradually increasing so that the average level of exogenous 
thyroid would build up to the desired value K’, would minimize the tran- 
sient effects and yield an observational means of determining the optimum 
rate of thyroid intake. The former method would certainly produce stabili- 
zation in a shorter time, and has been described by Gjessing (Joc. cit.). In 
this procedure, however, care must be taken to start treatment at an opti- 
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mum point on the oscillatory cycle and an extended period of preliminary 
study is required before treatment can begin. The second method allows 
treatment to be initiated at any time and avoids the danger of large peak 
values of thyroid concentration. 


The expenses of this investigation were defrayed, in part, by the Ada 
P. Kradwell Memorial Foundation and by a special gift of Mr. Arnold 
Frieder. 
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A general theory for the design of apparatus for the low temperature drying of frozen tissues 
has been developed. Using this theory it is possible to compute drying rates for different mate- 
rials in different apparatuses and to design an apparatus suitable for the required purpose. It 
is shown that for the drying of animal tissues in the laboratory, where the area is small, the re- 
quirements of the apparatus are not particularly critical. The same theory applied to the dry- 
ing of biologicals and plasma indicates that the design is exceedingly critical and suggests direc- 
tions in which existing apparatus might be modified. 


Introduction. In a preceding paper the general theory of the vacuum 
drying of frozen tissues was developed (Stephenson, 1953). In this paper 
the theory will be applied to the design of drying apparatus. 

Efficiency of drying apparatus.t We will define the efficiency of a drying 
apparatus as the ratio of the minimum possible drying time to the actual 
drying time. We have the general relation 

dx _dm_K 


Pee eT) aa 


SPs 
dt dt a 


(2eRT,) 7?" 


where 


w is the water content of the frozen tissue or solid, 

dx/dt is the rate of recession of the interface separating the dried and undried 
portions, 

dm/dt is the rate of evaporation per unit area, 

K is the general coefficient of net heat transfer to the evaporating surface, 
per unit area, 

L is the latent heat of evaporation of ice, 

T. is the average absolute temperature of the surrounding environment, 

T, is the absolute temperature of the subliming interface, 


* Aided by a grant from the Wallace C. and Clara A. Abbott Memorial Fund of the Uni- 
versity of Chicago and from the Commonwealth Fund of New York. 

+ Equations (1), (2), and (3) are derived in the preceding paper (Stephenson, Joc. cit.). The 
factor a has been omitted because it is dimensionless and equals unity. 
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P, is the saturated vapor pressure of ice at temperature 7, 

R is the gas constant per gram of water vapor, and 

fis the probability that an evaporated water molecule is removed from the 
system before it recondenses on the interface. All units are in the c.g.s. system. 


Also 
=—+—-l, (2) 


where 


fs is the probability that an evaporated water molecule reaches the exterior 
surface of the tissue before it is recaptured, and 

fe is the probability that a water molecule which reaches the exterior surface 
is removed from the system before it re-enters the dried portion of the tissue. 


From (1) we have for the time required to dry a thickness x; 


ne Ds x *t (2aRT,) 1/2 
swf Say —T) =wf ack Geena (3) 


As has been shown in the preceding paper, f., f-, and K/L can be deter- 
mined as functions of the thickness of the dry shell and (3) can be inte- 
grated. In order to obtain the minimum drying time, we assume that the 
temperature drop is negligible and that f. equals one. Then 


_ (2eRT,.) 2 f7* dx 
bin —- P. —f ig ’ 


where P, is the saturated vapor pressure at temperature T,. 
According to our definition, in order for a drying apparatus to have a 
good efficiency two conditions must be simultaneously satisfied 


T.=T. (5) 


(4) 


and 


—«1 or Jeae 1. (6) 


Under some conditions it may be impossible to satisfy (5) and hence theo- 
retically it would be better to define fin as the drying time when fe equals 
one and K is a theoretical maximum. However, this definition would be 
useless unless the theoretical maximum of K could be determined. 

When (5) and (6) are satisfied we have 


fey 


efficiency Y 1 — 7.27, (7) 


and, since 
Prfs=Pefe=Pif (8) 


APPARATUS FOR DRYING FROZEN TISSUES 25 


where P, is the vapor pressure at the exterior surface of the dry shell, 


nits da 
efficiency = 1 ry (9) 


8 


In general the efficiency of a drying apparatus must be measured either 
by determining f, and f, experimentally (Stephenson, Joc. cit.) or by meas- 
uring P, and 7, and so determining P./P,. In some important cases, how- 
ever, f. can be computed, and, using known values of f,, the efficiency can 
be computed. 

The reader cannot fail to notice that the efficiency of a given drying ap- 
paratus is always relative to the particular drying situation and in the 
abstract is non-existent. In some ways it would be best to abandon the en- 
tire concept of efficiency. However, since it seems firmly implanted in the 
literature, we have tried to define it in a useful way. It should never be 
forgotten that a large change in either f or K alone may have remarkably 
little effect on the drying time and so on the efficiency. In order to shorten 
the drying time appreciably the resistance to vapor flow 1/f must be de- 
creased and also the heat transfer coefficient K must be increased. 

Computation of f,/f.: 1) Plane parallel evaporating and condensing 
surfaces. Here, if collisions with air molecules are negligible, f. equals one 
and f equals f.. The effect of air will be considered below. 2) Evaporating 
and condensing surfaces separated by long, straight tube of circular cross- 
section, free molecular flow.* Using the fact that total flow through the 
tube must equal total flow through the dry tissue shell, we have 


AP $07 A ag 
(2aRT,) 1? ce ks he Se 


where 


A is the area of the interface separating dried and undried parts of the tissue, 
a is the radius of the tube, 

1 is the length of the tube, 

P, is the vapor pressure maintained over the condenser, and 

T is the average temperature of water vapor in the tube. 


When f. <f., P; K P., and T, = T, combining equations (8) and (10) 


gives 
ee ee (11) 


* See E. H. Kennard (1938, chap. VIII) for discussion of free molecular and diffusive gas 
flow in tubes. 
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3) Evaporating and condensing surfaces separated by a long tube of cir- 
cular cross-section, diffusive flow. Here we have 


AWE if = aa 2 2 12 


where 7 is the coefficient of viscosity. Under the same assumptions as for 
the derivation of (11), we obtain 


lice 16% /RT VATA 
7=| 7 21 a | ; Cy) 


In order to decide whether (11) or (13) is applicable, we must compute 
the mean free path in water vapor at pressure P, and temperature 7, and 
compare it with the radius a. If we take the mean free path \ at atmos- 
pheric pressure and temperature 150° C. to be 4.2 X 10-* cm., we have 
for the mean free path in cm. (Kennard, 1938, p. 149) 


7 beige 
h=6.5%10-2——. (14) 
P. 
or, if we ignore the variation of mean free path with temperature, 
PALS F 
— 3 
A=3 X10 FP. (15) 


Pressures in (14) and (15) are in mm. of Hg. For given values of efficiency 
f./f, temperature 7’, and load factor (Af,/), (11) and (13) can be used to 
compute the necessary radius. Curves computed in this way are shown in 
Figure 1. 

Effect of air pressure on efficiency. If the partial pressure of air is suffi- 
ciently low it will have negligible effect on the rate of drying. For example, 
at a pressure of 10-4 mm. of Hg the mean free path between collisions of a 
water molecule with air molecules is about 50 cm. Obviously, collisions of 
a water molecule with air molecules are of negligible frequency in com- 
parison with collisions with the dry shell, the walls of the drying chamber, 
and with other water molecules. In contrast, at atmospheric pressure the 
distance between collisions of a water molecule with air molecules is about 
5 X 10° cm. and collisions of this type will be the most frequent. The 
problem is to determine the effect of air at partial pressures between these 
extremes. In particular, we wish to know at what pressure air between the 
exterior surface of the tissue and the vapor trap begins to reduce apprecia- 
bly the rate of drying. 

We have obtained an exact solution of the effect of air on the rate of 
drying only for vapor flow between plane parallel evaporating and con- 
densing surfaces. Here, if the temperatures of evaporator and condenser 
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remain constant, after an initial period of transient flow, a steady state 
will be reached in which the net flow of air molecules through unit area of 
any surface between and parallel with the evaporator and condenser will 
be zero, and in which the net flow of water molecules will have some con- 
stant value. Thus any volume element of the gas will have a net velocity 
toward the condenser. Superimposed on this velocity will be an additional 


— DIFFUSIVE FLOW AvTe-10°C ¢ 
$2 0nG 2 =01 
5 Ola 
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E, -30° C8 =001 
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Ficure 1. Relation between efficiency, radius of tube connecting drying chamber and vapor 
trap, temperature and load factor f,Al. Here efficiency equals 1 — fe/f-, according to equation 


(7). 


diffusive flow. Therefore we have as differential equations describing the 
flow through any surface per unit area 


f dn; 


—D——+im=0, (16) 
dx 
Oe ar He 17) 
hee prt Sasa, 
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where 
D is the coefficient of diffusion of water vapor and air in cm.?/sec., 
nm, is the concentration of air molecules, 
mo is the concentration of water molecules, 
d is the velocity of flow, 
F is the mass rate of flow per cm.’, and 
ms, is the mass of a water molecule.* 


Since partial pressures are proportional to concentrations, the sum of 
the concentration or pressure is the same everywhere, and 3/D is inde- 
pendent of total pressure, we may solve equations (16) and (17), 
obtaining a eh 

d(x— ; 
Be PSE oe a8 (18) 
where 

p is the partial pressure of air at a distance « from the evaporating surface, 

p, is the partial pressure of air at the condensing surface, and 

1 is the separation of the condenser and evaporator. 


We relate (18) to the efficiency by imposing boundary conditions. We 
assume that the total pressure is everywhere the same. The pressure on the 
evaporating surface results from: 

1) The recoil momentum of those molecules emitted per unit area per 
unit time. This equals P,/2, where P, is the saturation pressure of water 
vapor at the temperature T, of the evaporator. 

2) The momentum absorbed from those molecules recondensing on the 
surface. If f, is the probability that a water molecule emitted from the 
evaporator reaches the condenser without being absorbed, this is 
(1 — f.)P,/2. 

3) The momentum absorbed from those water molecules emitted by the 
condenser which reach the evaporator. This is f,P,/2, where P, is the satu- 
rated vapor pressure at the temperature 7, of the condenser and f, is the 
probability that a water molecule emitted by the condenser reaches the 
evaporator and condenses. 

4) From air molecules. This is 

— oi 
ps = Pp, exp D (19) 


If we similarly analyze the total pressure on the condenser and assume 
Fm Tere (20) 


a for notation, equations (16) and (17) are identical with those used by A. J. Ede 
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we have, when we equate the pressures, 


eer ae Leeder. ‘ 
t= = ik = = at 8 r 
-AZt+s+sta= U-NZt+ S++ p. (21) 
This reduces to 
GQ—f) iP, = r— Pe= p.| 1—exp (™)]. (22) 
We have 
ae | tt 33 
—— Nt Ms P ’ ( ) 
where 
n, is the total number of molecules per unit volume, 
R is the gas constant per gram of water vapor, 
P is the total pressure, and 
T is the average temperature. 
Now, 


Pree an Sint ? poy tee Ble 
(2eRT.)/? (2aRT,)121 * (2aRT,) 7?" 


So we obtain, when we substitute (23) and (24) into (22), 


fl(P,—P,) RT? 
 (2aRT,) 2T'2 PD if (2, 


If we use the fact that RT*/PD is nearly constant, and take as physical 
constants 


F=j{ (24) 


(1-f) (P.—P,) = #,f1—exp [ 


D = .24 cm.?/sec. at 281° K. and one atmosphere, 
R = 4.62 X 10° erg./deg. per gram of water vapor, 
T = 281° K.., 

P = 1.01 X 10° dynes/cm.?, and 

1 mm. of Hg = 1.33 & 10% dynes/cm.?, 


we obtain 


P,— 1 — f) (P,—FP,) D IRN ee ee 
Pea exp [2.14 x 100 |, (26) 


where pressures are in mm. of Hg, temperatures in degrees absolute, 
length in cms. Equations (25) and (26) are similar to an equation devel- 
oped by P. C. Carman, although our method of derivation is different. 
Carman’s equation in a slightly modified form, which is identical with 
(25) and (26) except for notation, has been tested experimentally by H. 

* If T,~T,, the approximation is good and if T, > T,, P;/Ps< 1. In all these equations 


it should be remembered that vapor pressure is known as a function of temperature from data 
in vapor pressure tables. 
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Kramers and S. Stermerding (1951). That is, they have tested the 


equation 


Poa ae (27) 


where in our notation 
P =P, + p, = total pressure , 


Ilo. = P, , 
I’ = (1—f)P.+fP, , 
R/M=R, 
l=1, 
w= F,and 
Per; 


The subscript zero refers to 0° C. and one atmosphere. In general, agree- 
ment with experimental results was good. 

Equation (25) is easily transformed into a form which takes into ac- 
count the effect of the dry shell. From (25) we have, when P, K Ps, 


Pr a i f 
P, 1—exp(—kf)! fe 
where 
Pa 
oS dt ae 
b= 2.14X 10 ag. (29) 
If we write (2) in the form 
1 1 1 
Se es * 
PTS ee Se 
where 6 equals f/f,; write (8) in the form 
Patel. OF. 3 (32) 
and in (28) and (29) substitute P. for P, and f, for f, we obtain 
Pr Le te 
RNR ai SE it. EM EET eh lee hE 33) 
IF; PG f l ( 
g-—=* 1-— — 4-2 
- exp | — 2.14 10 Fgtraas fe| 
* When f;/fe & 1, we have 
fe 
6@~1—=—. 
fe OH 


Under these conditions, @ equals the efficiency as defined by equation (7). 
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Since from (30) 


i= 1 i aagac 
ea 
f 7, vi 7 1 Oli of Fi (1— 6), (34) 
(33) reduces to 
P,  1—exp(—k’6)’ oo) 
where 
/ RL 
hi 2. eae 
14 108 a (36) 
For ice, f, = 1 and @ = f.; hence (35) reduces to (28). 
| . & 
lo” 
(OF 
> | 
S) 10 
7 Ol 
u) 
a k=| 
s 
4 10 
0.01 
2 
10 
lo” 
4 
OO! lO 
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Ficure 2. Theoretical curves showing effect of air pressure on the drying rate of tissue and 
ice. The abscissae are values of #,/P; in equations (28), (35), and (55). The ordinate is f in (28), 
6 in (35), and 6’ in (55). It equals the efficiency, according to equations (31) and (7), under 
appropriate conditions. For (28) the parameter & is defined by (29), for (35) by (36), and 
for (55) by (56). 


It will be observed that the value of @ lies between zero and one, and 
serves as measure of the effect of air on the rate of drying. When #, equals 
zero, that is, when the residual air is negligible, 0, according to (35), must 
equal unity. On the other hand, when #,/P, > 1, k’0 is very small. In gen- 
eral, 9 must be computed as a function of the parameter k’ and p,/P.. 
Curves showing 0 as a function of p,/P, are shown in Figure 2, for values 
of the parameter &’ varying from 10‘ to 10~*. 
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Under certain circumstances useful approximate formulae for 6 can be 
obtained. When ’0 > 1, ie., R02 2.35, 


mc peae ss 37 
@=1— 5. (37) 
When 2’6 « 1, ie., 2’0 < 0.1, : 
= 38 
1S gts (38) 
P, 
and always 
— 2; 39 
6>1 P, (39) 
when k’@ < 1 
é> nate (40) 
i 


From equations (37) and (39) we can see that for 6 = 1, that is, for air 
to have negligible effect on the rate of drying, it is sufficient that p,/P, 
1. Certainly if p,/P; < .01, air will have negligible effect and very little 
if p,/P, < 0.1. It should be noted that so far this statement has been dem- 
onstrated only for plane parallel evaporating and condensing surfaces. 
The more general case will be considered below. As a glance at Figure 2 
or equation (38) will show, it is not necessary that ~,/P,; <1, for @ ~ 1, 
when &’ is sufficiently small. Ordinarily a small value of k’ implies a small 
value of f.. In other words, even a large relative value of the partial pres- 
sure of air has little effect on the drying rate provided its obstruction is 
small compared with the obstruction of the dry shell. 

It is sometimes interesting to know the conditions for diffusive flow be- 
tween the exterior surface of the dry shell and the condenser. This will 
occur when 


P, : 
po>iie, 210. (41) 
From (32) and (34) we have 
O fe 
P.=P, t= (1-0 (1 j,) 1P,> (1-8) Py. (42) 
Therefore 
Dirac ts Phdratie bate B 
FONE PT YE Ragk PTI ae 
So a necessary condition for diffusive flow is that 
1 
eres. (44) 
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or 
ROK1. (45) 


The argument as carried through does not establish (45) as a sufficient 
condition. However, when f, < 1, it is obvious that the approximations 
are mathematically good and that (45) is a sufficient condition. When fe 
approaches one, in ordinary drying situations k’@ <1 only when 6 <1, 
which again gives (45) as a sufficient condition. 

From (38) and (36) we see that when (45) is satisfied @ and 6f, are inde- 
pendent of P,, and depend only on f., #,, and J. (Temperature effects are 
of second order.) This is of course precisely what we would expect because 
(45) implies (41), as we have just demonstrated, which means that colli- 
sions of water molecules with each other are relatively infrequent relative 
to collisions of water molecules with the dry shell or air molecules. 

When @ <1, we obtain from (42) 


P.=>P, (46) 


and from (35) 


exp (— k’6) =1-3. (47) 


The general problem for two or three dimensions of the effect of air on 
flow in a drying system has not, so far as we know, been solved. An exact 
solution would depend on the solution of a differential equation analogous 
to (16), specifically, 

—div grad Dn,+div in, =0 (48) 
and 

—div grad Dn +div in: =0 (49) 
with imposition of appropriate boundary conditions. In any three-dimen- 
sional geometry, the solution of equations (48) and (49) is difficult, if not 
impossible. 

We suggest the following approximate solution. We will suppose that 
the resistance of the dry shell 1/f, and the apparatus without air 1/f, can 
be lumped into a single resistance 


hy Sawin (50) 


In order to compute the effect of the air we will introduce an additional 
hypothetical space between the drying apparatus and vapor trap to which 
the air will be confined, whose cross-section equals the actual average 
cross-section of the drying apparatus and whose length is the average dis- 
tance between the evaporating interface and the condenser. In this space 
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we will assume the flow obeys equations (16) and (17), that is, that its 
walls introduce no frictional effects. Under many circumstances these as- 
sumptions will be quite valid. Thus, when the ratio of air pressure to 
water vapor pressure becomes large, variations in pressure through the 
apparatus can be neglected and flow of water vapor will be analogous to 
flow of heat through a medium with constant heat conductivity and per- 
fectly insulating walls. 
Mathematically we will have 


Pe AGA 


see (E35) 


PaAafa= (Ps— Pa) A,f= 

where 

P, is the vapor pressure at the junction of the actual drying apparatus and 
the hypothetical added space, 

Ag is the cross-sectional area of the hypothetical space, which equals the 
average cross-sectional area of the actual apparatus, 

1/f. is the resistance of the hypothetical space to vapor flow, 

P, is the vapor pressure at the subliming interface, 

1/f is the combined resistance of the dry shell and the apparatus with no 
air present to vapor flow, 

A, is the area of the subliming interface between the frozen interior of the 
tissue and the dry shell, and 

6’ is defined by equation (54). 


Also 
ee Gee 2 
Po 1 et aia) 


(52) 


where & is defined by (29) with P, substituted for P., T for T., where T 
is the average absolute temperature of vapor in the apparatus, and / is 


equal to the distance between the tissue and the vapor trap. 
From (51) 


Agha dertemaelited he 
ir Bae Ee ME PS SEL 
and 
AAS 
aes Oe oe me tea) 


It should be noted that the definition of 6’ is such that when f. = 1, 
6’ = 1, and when f, = 0, 6’ = 0. 
Substituting (53) and (54) into (52) we obtain 


1— 9’ 


Texp(— 889’ Sie 


oe 
P 
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where 
Arf, 1 1 


k’’=9.0 X 104 : 
AgtA.f (2xRT)12 712? 


(56) 


where lengths are in centimeters, pressures in dynes/ cm.”, and tempera- 
ture in degrees absolute.* 
Equation (55) can be discussed in the same way as equation (35). In 
particular we have 
ees 
= P Git) 


8 


and, when k’’6’ < 1, 
1 


0’ = 5, (58) 


Af pr l 
Mee he eek th 
ig * Ae 5 eis (2ZakL) 3277 


Over the usual range of drying temperatures (0° C. to —50° C.) (58) be- 


comes 
' 1 
as (59) 


— Af ; 
1+ 100 p41 


where #, is in mm. of Hg and / is in centimeters. 

Effect of air on drying temperatures. When air is introduced into the dry- 
ing apparatus it will decrease the temperature drop of the subliming ice 
or tissue relative to the drying chamber. From (1) and (51) we obtain 


Tp Ba0rs Ve 
fg s ; ° 
te SK Aged (eRT VP eo 


In general for a given value of #,, (55) and (60) must be solved simulta- 
neously to give values of 7,, P;, and 6’. When (58) applies, 6’ is independ- 
ent of T, and, so, Ps. 

It will be observed that (60) can also be used to compute an experi- 
mental value of 6’, if JT, is measured as p, is varied. According to our 
theory, for a sufficiently small value of p,, 6’ approaches one and 7, be- 
comes constant. 

Curves for a particular apparatus (described previously, Stephenson, 
loc. cit.) showing T, as air pressure is varied are shown in Figure 3. The 
experiment was designed only to determine vacuum requirements for the 
particular apparatus and Macleod gauge. Hence, the data do not provide 
a critical test of the above theory. In particular, the experimental arrange- 

* In equation (51), P./(2RT,)*/? is the maximum sublimation rate of ice in gms./cm?. sec. 
(Fig. 4 in the preceding paper); 9.0 X 10‘ is the value of RT*/2/PD in c.g.s. units; (56) is trans- 


formed to the form of (29) and (36) by substituting the numerical values of R and of the factor 
for converting pressure in mm. of Hg into pressure in dynes/cm?. 
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ment was such that #,, the pressure at the condenser, was greater than the 
pressure measured on the Macleod gauge. Even without making this cor- 
rection, the data in general agree with equation (59). 

The effect of air on the rate of drying has been examined in such detail 
because in the biological literature, in particular, there appears to be con- 
siderable exaggeration of vacuum requirements (Malcolm, 1952; Packer 
and Scott, 1942; Scott and Hoerr, 1950). According to the above theory 
under no circumstances is there any necessity for the ratio of p,/P. to be 
less than .01, and frequently a much higher ratio can be tolerated. This 
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Ficure 3. Experimental curves showing the effect of varying air pressure on the sublima- 
tion temperature 7’, of tissue or ice. The average environmental temperature T, was 12°C. 
See previous paper for description of apparatus (Stephenson, Joc. cit.). 


fact has great practical significance because the requisite vacuum in a 
small system with a volume of a few liters can be obtained with a small 
mechanical pump. 

This theory does not imply that leaks can be neglected, because a leak 
not only lets air into the system but also a large additional amount of 
water vapor. 

A few examples will illustrate the application of the above theory. The 
minimum radius necessary for a certain efficiency* is a function of the 


* Here and in the following discussion efficiency means 1 — fe/fe in accordance with equa- 
tion (7), although the temperature drop T, — T, may not be negligible. 


—S st 
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temperature of the interface, the type of flow, and the composite “load 
factor” f,Al, where A is the area of the interface, 1/f, is the resistance of 
the dry shell, and / is the length of the tube [eqs. (10) through (13) and 
Fig. 1]. If we assume an interface area of 100 cm.2, a resistance of the dry 
shell of 1000, corresponding to a thickness of dry shell of .01 cm. for 
guinea pig liver, and a tube length of 100 cm.; we obtain a value of 10 cm.’ 
for f-Al. If the temperature of the interface is —30°C., the minimum 
radius required for an efficiency of 90 per cent (that is, f./f. = 0.1) is 2 
cm. if the flow through the tube is diffusive, and 2.3 cm. if it is free molecu- 
lar. From (15) we find 0.1 cm. for the mean free path. Hence the flow is 
diffusive. 

In order to determine the vacuum requirements for our apparatus, we 
use equations (55) through (59). From (56) a value of k”” of about 30 is 
found. Hence the assumption that k’’6’ <1 is not valid and required 
pressure is determined by (57). At —30° C. the vapor pressure of ice is 
about 0.3 mm. of Hg. This gives as the required value of air at the con- 
denser or vapor trap .03 to .003 mm. of Hg. 

It will be noted that this computation has been carried through for an 
interface temperature of —30° C. If this is the temperature of the drying 
chamber, the interface will be at some lower temperature as determined 
by equation (1). If we take K/L to be about 10-7 gms./cm”. deg. sec., 
which is in the range of an experimentally determined value for tissue in 
poor thermal contact with the drying chamber (Stephenson, Joc. cit.), we 
find the temperature of the interface to be about —42° C. This of course 
will affect the values for the radius, increasing it for diffusive flow to 2.5 
cm. for an efficiency of 90 per cent, and the vacuum requirement, chang- 
ing it to about 0.001 mm. of Hg. The mean free path will be about 0.4 cm. 
so the flow is still diffusive, although a transition to free molecular flow is 
beginning to occur. 

The value of 10 cm.* for the load factor f,Al is probably too high for the 
usual laboratory conditions. A value of 1.0 cm.* or less would probably be 
more typical and there would be a corresponding decrease in the required 
radius. For a given load and drying temperature the radius must be in- 
creased by a factor of 3.2 for diffusive flow and 2.2 for free molecular flow, 
in order to increase the efficiency from 90 per cent to 99 per cent. Thus, 
the additional increase in efficiency, which does not significantly shorten 
the drying time, is dearly purchased. On the other hand, if the radius is 
reduced very much below that needed for an efficiency of 90 per cent there 
is a significant increase in drying time. As long as the flow through the tube 
remains diffusive, the lower the drying temperature the larger the radius 
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needed for a given efficiency and load. For some loads, however, there is a 
transition to free molecular flow at low temperatures and the radius re- 
quired ‘for a given efficiency becomes independent of temperature. Thus, 
for a value of f,A/ of 1.0 cm.*, a drying temperature of —50° C., and an 
efficiency of 99 per cent, a radius of 2 cm. is required, assuming free molec- 
ular flow, and a radius of 6 cm., assuming diffusive flow. Since for this effi- 
ciency and temperature, the mean free path is about 10 cm., the assump- 
tion of free molecular flow is justified. If, as before, we assume a value of 
10-7 gms./cm.? deg. sec. for K/L and 10-* for f., we find that if the tem- 
perature of the drying chamber is —50° C., the temperature of the inter- 
face is —53.5° C., which has virtually no effect on the design requirement 
and very little effect on the rate of drying. 

As a final example we will compare the sublimation of ice and tissue for 
plane parallel evaporating and condensing surfaces, and then consider the 
effect of introducing an obstruction between the evaporator and condenser 
(Fig. 4). We assume that ice or tissue is in perfect thermal contact with 
a surface maintained at —30° C. by an appropriate supply of heat. Heat 
is conducted through the frozen part of the tissue or ice to the subliming 
interface where it is partly used to convert ice to water vapor and partly 
lost by radiation to the condenser. The water vapor then flows through the 
dry shell and intervening space to the exterior surface of the condenser 
where most of it condenses. The heat from the condensation and the radi- 
ation is then conducted through the ice on the condenser to a surface which 
is cooled to — 192° C. with liquid nitrogen. The basic equations for analyz- 
ing this drying situation are obtained by equating heat transfer through 
the ice layers, and the space between the evaporating and condensing 
surfaces. Thus we have 


out @e-T.) k(,-T,)__LgP,___—saL: SP, 
x — © wx  (2eRT,)1/2 (22RT,)17? Par 
+o(T3—T?).* 
Here 


Q is heat flow per unit area, 
xo is the total thickness of the tissue or ice, 
« is the thickness of the dried shell or sublimed layer, 


k is the coefficient of heat conductivity of frozen tissue or ice, which we have 
assumed to be .005 cal./cm. sec., 


* Equation (61) can be transformed into (1) by an appropriate definition of K and T, in (1), 
in terms of the variables appearing in (61). 
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Ficure 4. Comparison of sublimation of ice and drying of tissue for plane parallel evapo- 
rator and condenser and when the vapor flow between evaporator and condenser is obstructed 
by a long, narrow tube. In making computations, radiative heat transfer has been neglected 
and thermal contact of ice or tissue with the evaporating and condensing surfaces has been 


assumed ideal. 
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w is the water content of the tissue, 

T, is the temperature of the evaporator, —30° C. or 243° K. in our example, 
T, is the temperature of the subliming interface, 

T, is the surface temperature of ice on the condenser, 

T, is the temperature of the condenser, —192° K. or 61° K. in our example, 
L is the heat of vaporization of ice, which we assume to be 670 cals./gm., 
1/f is the resistance of the dry shell and any air in the intervening space, 

R is the gas constant per gram of water vapor, 
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Ficure 5. Curves showing variation of temperature of subliming ice surface T, as drying 
proceeds, for plane parallel evaporator and condenser. Uncorrected points refer to tempera- 


ture which would occur if heat transfer through ice on condenser was neglected and it was at 
absolute zero, 


P, is the vapor pressure at temperature 7.,, 
Ps is the vapor pressure at temperature T,, and 
ois the Stefan-Boltzmann constant (ice is assumed to radiate as a black-body). 


The vapor pressures are known as a function of the temperature from data 
in vapor pressure tables, hence equations (61) can be solved for T sand Ty. 
When these are known, the net rate of sublimation can be computed. 

In Figure 5 the temperature 7, is plotted for ice and guinea pig liver 
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for a total thickness of 0.5 cm. Radiative heat loss has been neglected. It 
is relatively negligible for ice and has very little effect on the drying rate 
of tissue, although for tissue it is much greater than the heat transfer by 
vapor flow. In the figure the uncorrected points refer to temperature T, 
which would occur if sublimation was into an absolute vacuum, that ist 
I’, was absolute zero and there was no ice layer on the condenser.When 
an obstructing tube is introduced, P, and the radiative term in (61) be- 
come negligible and a value of f is computed from (2) and (10) or (12).* 
Rates of sublimation for the obstructed and unobstructed flow are shown 
in Figure 6. 

Certain implications follow from the above examples. In the first place 
the obstruction has relatively little effect on the rate of drying of tissue 
but causes a considerable reduction in the rate of sublimation of ice. A 
substance such as penicillin or various other biologicals would be expected 
to dry very nearly like water, a substance such as plasma would be expect- 
ed to dry in an intermediate manner. Values of f, have not yet been deter- 
mined for such substances, although it is obviously important to do so. It 
is probable that considerable increase in the rate of drying of plasma could 
be effected by using a thin layer on the evaporating surface and molecular 
distillation apparatus. This is certainly true for biologicals. For a sub- 
stance which drys like ice, that is, which has a low solid content, the maxi- 
mum thickness can be computed once one establishes the permissible tem- 
perature drop relative to the maximum permissible drying temperature. 
Thus if the maximum permissible temperature is — 20° C. and one wishes 
to maintain drying at 90 per cent of the rate of maximum sublimation at 
—20°C., T, should be greater than —21°C., or the permissible drop 
(T, — T,) is 1° C. Using (56), the values of & and L for ice, f = 1, and as- 
suming P, is negligible, one finds for the maximum thickness of a, 7.3 X 
10-3 cm. For 50 per cent efficiency, a temperature drop of about 10 de- 
grees and a thickness of about .15 cm. is permissible. Obviously with such 
thin layers some sort of an intermittent cycle would be needed. However, 
a much greater yield of dried product per unit area of evaporator would 
result. 

The condenser temperature should be such that P, is a reasonably small 
fraction of P,, say, less than 0.1. Thus for a temperature of the subliming 
interface of —21° C., there is no benefit in having condenser temperature 
below —40°C. For T, equal to —30°C., T, should be equal to about 
—50° C. 

At —20° C. the vacuum requirements are somewhere between 10 and 
10-? mm. of Hg. 


* The approximations used in the derivation of (11) and (13) are not valid here. 
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Ficure 6, Curves showing sublimation rates of tissue and ice for plane parallel evaporator 
and condenser (unobstructed flow) and for evaporator and condenser separated by long tube 
as in Figure 4 (obstructed flow). As for Figure 5, uncorrected curve refers to sublimation rate 
which would occur if sublimation was into an absolute vacuum. The relation of the rate of 
sublimation to various limiting factors of vapor flow and heat transfer should be noted. If the 
total thickness of ice was slightly less, the final sublimation rate would be the maximum rate of 
sublimation of ice at —30° C, the evaporator temperature, 
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The meaning of ‘‘chloroplast reducing potential’”’ is examined in order to find out whether 
determinations of it might not allow one to circumscribe possible mechanisms of photosyn- 
thesis. One can define a physical situation where a measure of ‘‘chloroplast reducing poten- 
tial’’ could allow one to draw an important conclusion about the mechanism of photosynthesis. 
However, in order to achieve this, it is necessary to make very severe assumptions and to sup- 
pose that some very serious experimental difficulties can be overcome. 


Since Hill reactions (Hill, 1937, 1939; Hill and Scarisbrick, 1940a, 1940b) 
were discovered, the idea has been prevalent that it may be possible to 
derive a valuable conclusion about the mechanism of photosynthesis from 
a knowledge of the ‘“‘reducing potential” of illuminated chloroplasts. The 
nature of the Hill reaction and the considerable evidence (Clendenning 
and Ehrmantraut, 1950; Ehrmantraut and Rabinowitch, 1952; Hill, 1939; 
Holt and French, 1948; Macdowall, 1949; Spikes, 1952) which indicates 
that it is indeed a partial reaction of photosynthesis make this idea seem 
reasonable. However, the investigations on the ‘reducing potential” of 
illuminated chloroplasts (Gerretson, 1950; Macdowall, 1952; Spikes et al., 
1950) and the discussion relating to it which have to date been published 
have so far neither clearly defined the concept of “chloroplast reducing 
potential” nor shown how it could be of value in circumscribing possible 
mechanisms of photosynthesis. 

In the present paper an attempt is made to give an unambiguous defini- 
tion to a concept of “chloroplast reducing potential,” which, in principle, 
is experimentally determinable, and to show that it may be possible to 
draw an important conclusion about the mechanism of photosynthesis 
from a knowledge of it. It is necessary, however, to make rather severe 
assumptions and to suppose that some very serious practical difficulties 
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can be overcome in order to achieve this. The following analysis, therefore, 
probably has little value as a guide to fruitful experimentation, but it does 
demonstrate the great difficulties and ambiguities in the way of getting a 
useful measure of “chloroplast reducing potential.” 

A simplified, generalized scheme of photosynthesis. In order to be able to 
deal in more concrete terms with the possibilities of mechanism inherent 
in photosynthesis and the Hill reaction, the following generalized scheme 
(similar to that introduced explicitly by Hill, 1951a, 1951b, and certainly 
referred to implicitly by others) is introduced. The process of photosyn- 
thesis involves an over-all reaction: 


CO,-+H.0 — (CH20) +02, 


the energy for which is supplied by light. This reaction involves the re- 
duction of carbon dioxide by hydrogen coming originally from water. 
Carbon dioxide must be reduced to carbohydrate, and this reduction 
process can be expressed by the electrode reaction 


(CH:0) +H,0—CO,+ 4H*t + 4e E, = -42 


(Lipman, 1946). There is nothing implied in this about the mechanism of 
the reduction; it is rather a statement of its thermodynamic properties. 
Similarly, water must be oxidized to oxygen, which can be expressed as 


2H,0=—0.+ 4H++4e E;=—.815, 


with E’ referring to the electrode EMF (MacDougall, 1939) for standard 
conditions except that pH is 7. The two electrode reactions, if coupled, 
could go spontaneously only in the direction reverse that of photosyn- 
thesis. It is necessary, therefore, to postulate two other driving reactions: 
one to reduce carbon dioxide and the other to oxidize water. The scheme 
is represented in Figure 1, where energy directly or indirectly derived 
from light might be contributed either in the course of the transformation 
from QH to XH and X to Q, or directly to XH or Q, to drive the cycle 
of the photosynthetic mechanism. 

Since photosynthesis gives rise only to a net uptake of carbon dioxide 
and water and production of oxygen and organic matter, we may conclude 
that the detailed mechanism of photosynthesis involves cyclical use of 
the various “catalysts” concerned. In Figure 1 the generalized cycle of 
the hydrogen carrying catalysts is represented by the arrows surrounding 
the diagram. 

In Figure 2, where X and Y correspond respectively to X and XH, the 
large cycle of reactions represents the cycle of Figure 1 in more detail. 
Among the reactions entering into this cycle are reactions driven directly 
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or indirectly by light and the reaction which oxidizes water with the 
elimination of free oxygen. The R and RH are components of some redox 
couple upon which the reducing mechanism of photosynthesis operates; 
it may be the carbon dioxide-carbohydrate couple, the couple of reduced 
and oxidized Hill reagent, a redox indicator dye couple, an inert elec- 
trode, etc. 


—_—— X<-XH<———x_ E 34 
= 4 


(CH,0] +H,O CO,+4H*+4e\ E; 


: | 


xe 2H,0=0,+4H"+4e / £,=-815 
E % -815 


——— Q— QH ————> 
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Chloroplast reducing potential. The problem posed is one of determining 
the reducing capabilities of the couple X = Y in the intact biological 
system (chloroplasts) in a steady state. If one attempts to measure the 
reducing potential with either a redox indicator dye or an inert electrode 
immersed in the illuminated chloroplast suspension, then there are several 
possible results, as the following analysis demonstrates, and not all the 
possible results give useful information. 

For the sake of simplicity it must be supposed that the reactions around 
the cycle of Figure 2 are unimolecular. If catalysts are involved, this con- 
dition can only be fulfilled at low light intensities where, in photosynthesis 
and the Hill reaction, the rate is proportional to light intensity. Further- 
more, if one is to use the results obtained on a Hill reaction to draw con- 
clusions about photosynthesis, it is necessary that, except for the substi- 
tution of another oxidant for carbon dioxide, the Hill reaction mechanism 
is the same as the mechanism of photosynthesis, in regard both to reaction 
steps and to the rate constants of these steps. 

In addition to the reactions around the cycle of Figure 2 and the reac- 
tion of the X = Y couple with the system R = RH, there is the possi- 
bility that some other oxidizing couple, U = UH, foreign to the large 
cycle of Figure 2, can also react with R= RH; U = UH could, for in- 
stance, be part of the system involved in liberating free oxygen from oxi- 
dized water. 

The concentrations of components in the system of reactions of Figure 2 
are completely determined by NV + 2 steady state equations and three 
equations of conservation and definition. 

The steady state equations are 


Ry (Y) — k_y(X) = ky (RH) — k_, (RH) 
= kz (X) — k-x(A) =k, (A) —k-4(B) =... 
= ke (€) — kup (E+) 
= ky (M) — k-w (N) = kw (N) — k-n(Y), 


(1) 


where k_, and k, are abbreviations for k_y(UH)and ky(U) respectively, 
and k_y and ky are abbreviations for k_z(RH) and k,(R). 
The conservation and definition equations are 


So= (X) + CY) CA) +. 2+ CN) (2) 

Ro = (R) + (RH) (3) 
(RH) 

Tae % 


CHLOROPLAST REDUCING POTENTIAL 49 


The quantity Rp refers either to total amount of redox indicator dye added 
or to the surface area of the inert electrode exposed to the chloroplast 
suspension. 

The details of the computations involved in solving this system of equa- 
tions to give J as a function of the constants of the system are shown by 
Horwitz (1952). 

The final steady state equation relating I to the constants of the sys- 
tem is 


Ki Rol K2Ro 


SoK Rot SoKeRol = (boRol — k-wRe) | Ey + PE Ty 


+Ks], (5) 


where Ki, Ko, K3, Ks, and K; are complex constants. 

For the simple case where there are only three components, X, Y, and 
A, in the large cycle of Figure 2, one easily arrives at essentially the same 
result, except that the makeup of the K’s is simpler. By solving the three 
steady state equations of this simple system, one finds that 


_ [ku(RH) — k-.(R)] {k-vyRath-vk_-xthxka} 


2 ee ee a (6) 


x, 


and X and A are given by cyclic permutation of this result: Y ~ X —> 
A—Y. 
Making use of the conservation equation involving So, one finds 


[ku (RH) — k-u(R)] {Ki (RH) + K2(R) + Ks} (7) 
K,(R) + Ks (RH) 


So= 


where Kj, Ko, K3, K4, and K; are constants. By eliminating R and RH with 
equations (3) and (4), one easily obtains a result formally identical with 
equation (5). 

From this it is apparent that J, which is the quantity one would meas- 
ure experimentally in such a method of determining chloroplast reducing 
potential, is a function of Ro, the amount of redox indicator, and So, corre- 
sponding to the amount of chloroplasts used. Now, if, having found a 
redox indicator which comes to a partially reduced steady state with a 
suspension of illuminated chloroplasts, it turned out experimentally that 
the ratio J was a function of Ry and So, then there is little or nothing of 
consequence that one could conclude about the photosynthetic mecha- 
nism, except that it is capable of achieving the reducing potential repre- 
sented by the ratio 7, and possibly higher ones as well. 

However, if reactions like those whose rate constants are k, and k, 
and whose effect is to cause the cycle to undergo a net reaction in one di- 
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rection, are absent or negligible, then there will be a state of equilibrium 


and the steady state equation (5) reduces to 


Ky___ kykgka.. kw & 
K; k-yk—xk-a niece k_n 


I=- 


Thus in this case, I will be independent of the amount of redox indicator 
Ry added and chloroplast concentration So, and will depend only on the 
constants of the system. It is thus possible, in principle, to distinguish ex- 
perimentally between the ‘‘steady state” case, which can give us little in- 
formation, and the ‘“‘equilibrium”’ case. 

Practical difficulties, however, cause the clear cut distinctions enunci- 
ated in the last paragraph to lose much of their definiteness. For there may 
be experimental complications involved in determining that J is inde- 
pendent of Ry and So, since some of the k’s involved in ks and k; contain 
factors for light energy which might vary due to shading by redox indi- 
cator or chloroplasts themselves, and these variations in k’s would cause 
changes in J. In addition, since no estimate has been made of how much 
variation in J is to be expected from a given variation in Rp or So, it may 
be difficult to decide whether the equilibrium case obtains or J is changing 
very slowly with Ro or So. 

If an “equilibrium” situation is achieved, then it can be claimed that 
the reducing potential attained in this way is higher than the reducing 
potential obtaining in a natural photosynthetic steady state. This follows 
from the cyclic nature of the hydrogen transport mechanism represented 
in Figure 2. In a natural photosynthetic steady state there must be a nega- 
tive free energy change in going around the cycle in a clockwise direction 
from X to Y since there is a net reaction in this direction. In this manner 
of speaking, of course, one includes energy sources, like light energy, 
which drive the cycle, as reactants in the clockwise series of reactions 
from X to Y, whereas in an equilibrium this free energy change is zero, 
for there is no net reaction, which must involve an increase in the ratio 
(Y)/(X) to achieve the equilibrium. 

There are other possibilities of reaction schemes giving the “equilibri- 
um” criterion of an J independent of Ro and S». Thus, the &, and A, re- 
actions may be directly concerned with X and Y as shown in Figure 3, 
or the &, and k_, reactions might be concerned with any other step in the 
hydrogen transport cycle. In either case J is independent of Ro since, 
among equations determining these systems, (R) and (RH) occur only 
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in the equation expressing the steady state condition for the concentra- 
tions of R and RH 


Pe sero (Rn) er! oat) 
k—p(X) © 
That J is independent of S can be seen from the consideration that around 
the cycle of Figure 3 there are V + 1 steady state equations determining 
N + 1 ratios of concentrations of components around the cycle to concen- 
tration of, say, X. Introduction of the conservation equation involving So 


(9) 
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will not determine the ratios (and therefore not change J), but merely de- 
termine the absolute amount of each component around the cycle. 

Now the argument to demonstrate that even in these cases the ratio 
(Y)/(X) must be higher than in a natural photosynthetic steady state is 
not quite as rigorous as in the first case considered. An oxidizing couple, 
like U = UH, which provides a continuous supply of oxidizing power, 
could, in the absence of air, only be some intermediate between “oxidized 
water,” formed at one stage in the hydrogen transfer cycle, and free oxy- 
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gen. Presuming that the naturally occurring parts of the hydrogen transfer 
cycle of Figures 2 and 3 are removed from the cell undamaged with the 
isolated chloroplasts, then the constants k, and k_, should be no greater 
here than in a normally photosynthesizing cell, and the rate of breakdown 
of U to oxygen should be no smaller. Therefore, it is only by an increase 
in the ratio (Y)/(X) over that in normal photosynthesis (where there are 
two pathways available from Y to X: that by reduction of carbon dioxide 
and that by reduction of U) that a steady state can be maintained from 
Y to X only through reduction of U. One might be more satisfied of the 
validity of the supposition of an undamaged cycle if whole cells, such as 
have been used to give ordinary Hill reactions (Clendenning and Ehrman- 
traut, 1950; Ehrmantraut and Rabinowitch, 1950) were used, rather than 
isolated chloroplasts. 

On the carbon dioxide side of the cycle represented in Figure 1, there are 
two general possibilities of mechanism: 1) An X = XH couple with sufh- 
ciently high reducing potential to reduce carbon dioxide or its “exergoni- 
cally derived’’ fixation products; 2) Chemical or physical mechanisms 
acting on the components of the carbon dioxide-carbohydrate couple to 
lower its electrode EMF to a level where X = XH can reduce it. An ex- 
ample of such a mechanism is provided by S. Ruben (1943). 

Since the equilibrium case always gives a value for reducing potential 
higher than the natural one in steady state photosynthesis, one can con- 
clude, if the equilibrium reducing potential is less than sufficient to direct- 
ly reduce carbon dioxide, that one of the set of mechanisms under 2) in 
the preceding paragraph must be operative. Even this conclusion, how- 
ever, must be based on an assumption about the photosynthetic mecha- 
nism. For, the principle that, in a complex set of reactions resulting in the 
net addition of four hydrogens to carbon dioxide, the reducing agent 
causing this reduction must have a higher reducing potential than that 
of the carbon dioxide-carbohydrate couple is only valid if the net of four 
hydrogens added actually represents the only hydrogens transferred. 
Thus, it is not inconceivable that during the transformation from carbon 
dioxide to carbohydrate, say, eight hydrogens are added, four of which 
are removed in reactions which are not simply the reverse of the addi- 
tion, but result in a net storage of energy. 


LITERATURE 


Clendenning, K. A. and H. C. Ehrmantraut. 1950. ‘‘Photosynthesis and Hill Reactions by 
Whole Chlorella Cells in Continuous and Flashing Light.” Arch. Biochem., 29, 387-403. 
Ehrmantraut, H. C. and E. Rabinowitch. 1952. ‘Kinetics of Hill Reaction.”? Arch. Biochem. 

Biophys., 38, 67-84. 


CHLOROPLAST REDUCING POTENTIAL eS, 


Gerretson, F. C. 1950. ‘Manganese in Relation to Photosynthesis II. Redox Potentials of 
Iuminated Crude Chloroplast Suspensions.”’ Plant and Soil, I, 159-93. 

Hill, R. 1937. “Oxygen Evolved by Isolated Chloroplasts.’’ Nature, 139, 881-82. 

. 1939. “Oxygen Produced by Isolated Chloropiasts.’’ Proc. Roy. Soc. Lond., B 127, 

192-210. 


. 1951a. ‘‘Oxidoreduction in Chloroplasts.” Adv. in Enz., XII, 1-39. New York: Inter- 

science. 

. 1951b. “‘Reduction by Chloroplasts.’’ Symp. Soc. Exp. Biol., V, 222-31. Cambridge: 

Cambridge University Press. 

and R. Scarisbrick. 1940a. ‘“‘Production of Oxygen by Illuminated Chloroplasts.’’ Na- 

ture, 146, 61-62. 

. 1940b. “The Reduction of Ferric Oxalate by Isolated Chloroplasts.”’ Proc. 
Roy. Soc. Lond., B 129, 238-55. 

Holt, A. S. and C. S. French. 1948. ‘Isotopic Analysis of the Oxygen Evolved by Illuminated 
Chloroplasts in Normal Water and in Water enriched with O18.”’ Arch. Biochem., 19, 429-35. 

Horwitz, L. 1952. ‘‘Investigation on Hill Reactions of Isolated Chloroplasts and Whole Cells 
and Their Relation to Photosynthesis.’’ Thesis. University of Minnesota. 

Lipman, F. 1946. ‘‘Metabolic Process Patterns.’’ Currents in Biochemical Research, pp. 137-48. 
(Ed. D. E. Green.) New York: Interscience. 

MacDougall, F. H. 1939. Thermodynamics and Chemistry. New York: John Wiley and Sons. 

Macdowall, F. D. H. 1949. ‘‘The Effects of Some Inhibitors of Photosynthesis upon the Photo- 
chemical Reduction of a Dye by Isolated Chloroplasts.’ Pl. Physiol., 24, 462-80. 

. 1952. ‘“The Reducing Potential of Illuminated Chloroplasts.’’ Science, 116, 398-99. 

Ruben, S. 1943. ‘“Photosynthesis and Phosphorylation.’’ Jour. Amer. Chem. Soc., 65, 279-82. 

Spikes J. D. 1952. ‘Stoichiometry of the Photolysis of Water by Illuminated Chloroplast 
Fragments.”’ Arch. Biochem. Biophys., 35, 101-9. 

, R. Lumry, H. Eyring, and R. E. Wayrynen. 1950. ‘‘Potential Changes in Suspensions 

of Chloroplasts on Illumination.”’ Arch. Biochem., 28, 48-67. 


RECEIVED 4-10-53 


BULLETIN OF 
MATHEMATICAL BIOPHYSICS 
VOLUME 16, 1954 


TRANSIENT PHENOMENA IN CAPILLARY EXCHANGE 


H. D. LANDAHL 


COMMITTEE ON MATHEMATICAL BIOLOGY 
THE UNIVERSITY OF CHICAGO 


The case of a vessel, which supplies a region through which it passes with some substance, 
is considered for the situation in which the permeability is the limiting factor. Diffusion parallel 
to the vessel is neglected. The substance may, however, be consumed proportional to its con- 
centration in the inner or outer region. A solution is given for the case in which the input is an 
arbitrary function of time. It is suggested that the results may be applied in some cases to data 
on the injection of substances into blood vessels, or they may be applied to the transient effects 
in the case of vapors or gasés passing through the respiratory passages. 


The problem of exchange in capillaries has been treated recently by 
G. Schmidt (1952). A special case of this problem has been treated by 
W. C. Sangren and C. W. Sheppard (1953). In this paper only the per- 
meability of the capillary wall was taken into account and no reaction was 
considered to take place in the tissues supplied by the capillaries. The 
initial conditions were not easy to meet since they required rapid injection 
without, however, altering the velocity. It is the purpose here to consider a 
somewhat more general case than considered by Sangren and Sheppard, 
the boundary conditions assumed being closer to the conditions that are 
likely to be the case in an experimental situation. It is not necessary to 
know the initial distribution along the length of the capillary. 

Consider a system consisting of a vessel and its surrounding tissue in 
which the permeability of the wall is the limiting factor for the substance 
S being considered, delays due to diffusion in the surrounding tissue being 
neglected. The length of a capillary is generally much greater than the 
“diameter’’ of the region which is supplied by it. We shall therefore neg- 
lect any transport parallel to the vessel in the region about the vessel. Let 
the velocity of flow through the vessel be a constant. Let r be the radius 
of the vessel, v, the volume associated with each unit length of the vessel, 
h and k the permeability coefficients for inward and outward flow, F the 
concentration of S within the vessel, and G the concentration of S in the 
surrounding region. We shall suppose that the substance may be con- 
sumed proportional to its own concentration in the surrounding region, 
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the coefficient of proportionality being a. We shall consider the situation 
in which initially F and G are both zero. Let the input at x = 0 be a 
given function C(t) = F(0, t). We shall suppose that the injection produc- 
ing C(t) does not appreciably alter the velocity V. The differential equa- 
tions and boundary conditions corresponding to this problem are then as 


follows: 

86 = bAP— AG ~aG, (1) 
OF OF 
Ahsan) ae es 2 
a bBF + BG Fea (2) 

F(x,0) =G(«,0) =0, _F(0,#) =C(), (3) 
ees ene ed (4) 

Vs k r 


We may next introduce 
i (#8) =8(F (x, )}= fo eR (a, dt, 
0 


and similarly define g(x, s) and c(s) as the transforms of G(x, ¢) and C(¢). 
Since in a physical problem the real part of s will be positive, we find from 


(1) and (2) 
sg=bAf—(A+a)g, (5) 


sf= —bB/+Bg—v 9f. (6) 


Solving for g in (5) and substituting for A in (6) we find on rearranging 


a janie ite) Gea: 
ax (0B+ s stA+ta/V' eed 
so that if 
bAB 
a=, (8) 
f(a, s) = eR M gH) EA (ot Ad a)e(s). (9) 
st+tA+a ‘ 
Hence on inverting we find 
(a/s+ A+a) 
F(a, 1) =e—(@B2/V)| @-1) SO CQ 
(x, i) =e [2 a} g [(stAta)e(s)}] (10) 


in which 2" is the inverse Laplace transform and the asterisk indicates a 
convolution integral. If we wish C(0) to approach a finite value Co, we 
may let C(t) = Cot/e for ¢ < ¢ and C = C(#) for ¢ > € and then let € be- 
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come very small. This introduces an additional term, which, when added 


to (10) after the latter has been written explicitly, becomes 


wt Ux; Tt) = Coe e—tatoc-tl) Tol a(t =) ] 


$e-omin) [OT A +a)C (1-4-7) +0, (:-5-7)] (11) 
X elt") (2 Var) dr, 


in which C; denotes the derivative of C with respect to time and Jy denotes 
the zeroth imaginary Bessel function. If one integrates the above expres- 
sion by parts, the result is 


F (x, t) =c (1-4 oye (6B2/V) 4. e— epalv) fo ask. c (1-4-7) 


X e—(Ata)rI, (2 Var) dr, 


in which J, is the derivative of Zo with respect to r. 

Ifa=0,6=1, B=p, A = p’, and if C(é) is zero everywhere except 
at ¢ = 0, but its integral is M, then the above is the same as equation (7) 
of Sangren and Sheppard (loc. cit.). 

We consider in somewhat greater detail the case in which C(#) is a 
constant C. Then, for not too large values of the time after starting, we 
find on expanding and introducing ? = ¢— x/V and A +a=y 


F(x, 0) =Coe02/)| Ty (2 al) Lael — yt— e-11) 


73 RP (13) 
a 
2 N+ (a 3y) gt (a? 6ayt 12) oat t---]- 
If t is very large, then it is easily shown that 
F(x, ©) =Coe(B2/(440)¥] , (14) 


If the substance S is consumed proportional to its concentration within 
the vessel at a rate GF, a being the coefficient of proportionality, then the 
above expressions also represent the solution for this case if the quantity 
dB, occurring explicitly, is replaced by bB + 4. 

If the vessel length is ZL and we wish the rate at which the substance S 
is leaving the vessel, this value 17?V F(L, ¢) can be found by setting « = L 
in the above expressions (11) through (14) and multiplying the result by 
arV. 
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It may be pointed out that the coefficient of consumption a may be due 
to disappearance throughout the volume or it may be that we consider 
the vessel to be a vessel carrying air, the vessel wall now being the sur- 
rounding medium and at some average distance the blood vessels remove 
the substance proportional to the concentration. In this case these results 
may be useful in the consideration of the removal of some vapors or 
gases from the respiratory tract (cf. Landahl and Herrmann, 1950). 

In some situations it may be convenient to saturate the system with the 
substance S, then at ¢ = 0 the system is perfused with blood not con- 
taining S. The solution in this case is just the complement of the case 
above in which C is a constant and a = 0. This can be seen if it is noted 
that initially G(w, 0) = 6 F(a, 0) = bCo, Co being a constant. Then intro- 
ducing bC) — G(x, #) and Cy — F(a, ¢) as new variables one finds that 
these variables satisfy (1) through (3) if substituted for G and F. 

If the substance S is consumed within the vessel at a rate proportional 
to the concentration, the coefficient being &, the solution for this case is 
given by (11), (12), and (14) in which 6B + & replaces 6B wherever the 
latter occurs explicitly. 


The author is indebted to Dr. G. W. Schmidt for reading and discussing 
the manuscript. 
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An approach to quantitative work on optimal systems is considered. Desired optimal prin- 
ciples are utilized in constructing a hypothetical system similar to the organ system con- 
sidered. A comparison of this constructed system with the anatomical system then gives an 
indication of the importance of the optimal principles in the form or function of the organ sys- 
tem considered. 

These ideas are applied to the mammalian vascular system, and limiting values are obtained 
for some of its important component parts. The constructed system gives good agreement with 
anatomical data for vessel radii, lengths, and hydrodynamic resistance to flow. 


Introduction. There have been several mathematical considerations of 
organic form in biological literature. Foremost among these is D’Arcy 
Thompson’s monumental treatise On Growth and Form (1917). He ana- 
lyzes many individual forms and extends his mathematical treatment to 
homeomorphisms between closely related forms in his now famous theory 
of transformations. Thus we find his illustrations graphing the profiles of 
the porcupine fish Diodon into the sunfish Orthagoriscus (Thompson, Joc. 
cit., p. 751) and the skull of Eohippus into Equus (Thompson, Joc. cit., 
p. 765). 

This method of applying mathematics to morphology has many pit- 
falls. Waddington (Zuckerman, 1950, p. 512) has discussed the application 
of topology to studies on form: 


From the strictly geometrical point of view, each new isolated cavity appear- 
ing within a body, such as a blood island in an embryo, or a food vacuole in an 
amoeba, raises the degree of complexity of the form by one unit, whereas a 
new external excrescence, even if it is as functionally important as a limb, leaves 
the degree of complexity unchanged. It is obvious in this case that we cannot 
pass simply from the mathematical relations to statements about biological rela- 
tions. We should require for that a specifically biological analogue of topology— 
a new branch of science in which relations of the logical kind used in normal 
topology were set up between entities whose definition was essentially biologi- 
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While opening new regions for the application of mathematics to biol- 
ogy, D’Arcy Thompson was led astray through a natural desire to apply 
in toto the methodology and concepts of an abstract science to a very con- 
crete science. Form can be analyzed mathematically, but are we to con- 
sider an inclusion of a different order of magnitude from a limb? We have 
not yet considered the relative functional utility of the limb and the in- 
clusion. Thus we must temper our mathematics with practical considera- 
tions. Rather than abstract mathematics let us attempt an application of 
engineering principles. 

Constructive work in this direction has been done by N. Rashevsky 
(1948). In the last few chapters of his book he considers problems of gen- 
eral organic form. Of particular interest to us are his considerations on 
tree and plant form and the appendages of animals. All these topics are 
considered from the viewpoint of functional utility. 

As our contribution to considerations of form the following is suggested 
as a fertile approach to the problems involved. 

1. A careful survey of present knowledge concerning the anatomy and 
functions of the organism or organ system to be studied. 

2. By considering the functions of the organism or organ system we 
must design a practical system or organism to accomplish these functions; 
call this the unmodified optimal system. 

3. We must compare our designed system with the actual system in 
existence. 

4. In so far as the designed system differs from the actual system we 
must consider modifications of our principles of design which will make the 
designed system more closely resemble the actual system. 

5. We must justify these modifications. 

6. We must investigate the possibility of other optimal systems; and, 
what is less probable, the possibility that our optimal system could have 
been derived from different principles of design. 

It is not sufficient that we criticize the classical approach to morphologi- 
cal study and substitute our own methodology. We must now justify the 
methods we have chosen. This can best be done by showing how they will 
increase our knowledge of the fundamental processes of living organisms. 

The form of the mature organism or organ system is directly dependent 
on and is a result of early embryological growth. This growth is deter- 
mined partly genetically and partly environmentally. What we mean by 
genetic determiners of growth is the existence of certain physiological 
mechanisms influencing growth that are transmitted from generation to 
generation. A specific example of such a mechanism is the genetically de- 
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termined production by the cell of a specific chemical substance that will 
influence growth. Such a substance might be an enzyme, or, as a more 
complex example, an enzyme system utilized in some specific way in the 
metabolic pathways of the cell. With such an enzyme system present 
utilizing some specific substance for metabolism, we further postulate a 
source of this substance outside of the cell. If the reaction this enzyme or 
enzyme system affects were to be necessary for cellular growth and/or 
division we would find no such growth were this extracellular substance 
not present. In general, there would also exist a gradient of growth such 
that the greatest growth would take place in regions of highest concentra- 
tions of this substance. Let us now postulate another enzyme system nec- 
essary for growth utilizing another extracellular substance. If growth de- 
pends on both systems being operative, growth will be along the gradient 
of some composite function of the two concentrations. 

Let us now consider a very complex field consisting of many gradients 
with growth determined by some complex combination of these gradients. 
An alteration in one of the enzyme systems responsive to one or several of 
these many gradients would then influence growth to a minor extent. This 
we might imagine as the process determining minor genetically deter- 
mined variations in growth within a species. Naturally, there will also exist 
environmentally determined variations which are of no concern to us here. 

This is a very simplified picture of growth determiners and growth 
gradients. In more general terms growth can be along observable chemical 
or physical gradients or what we may call random growth along no 
gradients we can subject to analysis. 

We may now proceed to consider a more complex subject: evolution and 
stable forms. To clarify our thinking let us consider as a specific example 
the form of the aorta, a part of the vascular system which we will later 
consider in more detail. Let us start with an organism with a given mass 
and weight, etc., and an aorta of given radius. Hydrodynamic considera- 
tions show that the pressure drop across the aorta is negligible compared 
to the pressure head necessary to force blood through the vascular system, 
and the velocity of blood in the aorta is such that the Reynolds’ number 
for turbulence is far below the critical value. During the course of genera- 
tions this same animal evolves into a form with an aorta of smaller radius. 
We find, as before, a non-critical Reynolds’ number and low resistance. 
The animal will be more efficient (and so better able to survive) for a) it 
will have less vascular tissue to maintain and b) with the increased cardiac 
output of a stress situation its venous return will be faster, a consequence 
of the elastic properties of the smaller aorta. Let the aortic radius keep 
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decreasing in magnitude and the Reynolds’ number will eventually be- 
come critical. The resistance will increase and a point will be reached at 
which the economy of maintaining a smaller aorta will be more than offset 
by the increased cardiac work necessitated by the high aortic resistance. 
By considering these factors we see that an optimal aortic form will evolve. 
This form will be optimal in the sense that, considered as part of a hydro- 
dynamic system, it means the least energy dissipation for the organism. 
The growth determiners determining aortic form in the individual must 
therefore be responsive, in the long term sense of evolution, to the hydro- 
dynamic factors discussed above. 

Ontogenetic development with an optimal organ system is a more 
subtle problem. The final form of the individual is present in latent form 
in the genetic makeup of the single cell from which development begins. 
Referring to the example used above, in so far as the genes responsible for 
the aortic form will give rise to an optimal form, we can say that they 
themselves are an optimal system. Considering the genes at the egg stage 
without reference to their later determining properties, it is difficult to say 
to what extent they then are an optimal system. However, we need not 
concern ourselves with this problem at present. 

The growth of the individual is no more than the sum of the growth of 
its parts. Thus, the growth gradients previously discussed will, when their 
total effect is considered, determine the growth of the entire organism. In 
a general sense we also see in what manner an optimal form of an organ 
system may be achieved through the random variations of evolution. 

There is a formal resemblance between the above optimal growth 
gradients and the modified optimal principles of design we have men- 
tioned in points 2. and 4. of our method. This resemblance lies in the fact 
that our constructed pencil and paper model is a direct result of the prin- 
ciples we have chosen in the same sense that the organ structure is a 
direct result of the growth gradients of the organism. The fertility of our 
method or approach centers around this similarity. 

The criticism of the classical approach to morphology outlined at the 
beginning of this paper leads to the question: If this approach is sterile, 
what more dynamic method of study would we offer as a substitute? We 
have discussed a substitute method. The reason for its hoped for success is 
that from our principles of design we can now arrive at an idea of the 
growth factors influencing growth of the system. Thus from pure consider- 
ations on form in relation to function we find a tie-up of morphology with 
the studies of growth and embryology. In other words, a true study of 
growth and form both as separate and as related disciplines. 
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The true test of any method of study lies in its practicality. As a first 
application of our method we will study the structure of the mammalian 
vascular system. 

The vascular system. Unicellular and simple metazoan forms can be sup- 
plied with food through simple diffusion. With the development of higher 
forms and specialization of lungs and food acquiring and storing organs 
came the development of the vascular system. This system functions in 
distributing metabolites to and draining waste products from the various 
tissues. Highest on the evolutionary scale stands the mammalian vascular 
system, for which we have a large store of anatomical and physiological 
data. 

The mammals cover a large range of sizes; from the mouse to the man 
to the whale, a mass variation of a million fold. The surface area of an 
animal is a power of its mass, and metabolic rate depends on surface area. 
Thus, as M. Rubner (1883) has shown, the mass determines metabolic 
rate; and, in fact, the metabolic rate varies approximately as the two-thirds 
power of the mass. The metabolic rate determines the oxygen needs of the 
tissues, and oxygen is transported to the tissues via the blood. The oxygen 
carrying capacity, and, in fact, the general physical characteristics of 
blood do not vary much among the mammals. For example, the red blood 
cells of the mouse, man, the elephant, and the whale are very similar in 
size. The mass of the animal determines the metabolism, the metabolism 
determines the oxygen need, and the oxygen need determines the blood 
flow per second (cardiac output), since all mammalian blood is about the 
same in physical characteristics (i.e., oxygen carrying capacity). Thus, 
from very general data and the known oxygen carrying capacity of the 
blood, we can relate cardiac output to animal size. Let us now consider 
other aspects of the vascular system. 

The various parts of the mammalian vascular system are the heart, the 
arterial system, the capillaries, the venous system, and the blood and 
blood forming organs. The arterial system and capillaries will be our sole 
concern. The venous system returns the blood to the heart, but does not 
play any great pressure or flow regulating role since it is a system of low 
resistance, thus we will not be concerned with it. Similarly the blood form- 
ing organs are of no concern to us. The physical characteristics of the 
blood we will accept as given for all the mammals. The heart is of great 
importance in the dynamics of the circulation but we will not treat it here. 
We will consider the arterial system as a passive system of vessels, with a 
pressure head maintained by the heart at one end. 

In treating the arterial system and capillary bed we will first investi- 
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gate the optimal size of the largest and smallest vessels and then will try 
to design an adequate system of connection. Thus our problem has been 


separated into the following parts. 


I. Size of the aorta. 
II. Capillary dimensions and the volume supplied by one capillary. 
III. The connecting system between the aorta and the capillaries: 
A. The structure of the branching system. 
B. Relative diameters of vessels at a branching. 
C. Relative lengths of the branches. 
IV. Total resistance of the system to flow. 
V. Comparison of our constructed system with an actual anatomical 
system. 


I. Size of the aorta. We have given a brief introduction to the nature 
of this problem in the introduction. Since the aortic resistance to blood 
flow is negligible compared to the total resistance of the system, we want 
to set a lower limit on its radius such that turbulence will not develop. 
The Reynolds’ number, call it Re for blood flow, has been found to be 
about 1100 for the beginning of turbulent flow (Coulter and Pappen- 
heimer, 1949). 

Let C be the blood flow in cm.’ per second (or cardiac output). Then, if 
ra is the radius of the aorta and @ the average velocity of flow in the aorta, 


we find 
G 


v= 2 
2 
rr 


(1) 


The Reynolds’ number for viscous flow in a smooth tube such as the aorta 
is defined by the equation 
Re =", (2) 


where 4 is defined as above, 7 is the viscosity of the fluid, o is the density of 
the fluid, and R is the radius of the vessel. Solving for the radius we obtain 


nRe 


ov 


Re 


(3) 


Substituting (1) in (3) and inserting the value 1100 for Re we obtain as our 
equation limiting aortic size 


eee 
met thea (4) 
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IT. Capillary dimensions and the volume supplied by one capillary. The 
main function of the capillaries is the exchange of metabolites between the 
blood and the tissues. This exchange occurs across the capillary mem- 
brane. Thus we want the surface area for this exchange to be as large as 
possible. This can best be accomplished by having a small capillary radius 
and many capillaries. Since we assume the blood characteristics as given, 
the lower limit to capillary radius is the size of the red blood cell. The 
blood should spend enough time in the capillary so that adequate metabo- 
lite transfer occurs. Since oxygen diffuses very quickly across the capillary 
membrane we are limited by the time it takes a slowly diffusing metabo- 
lite to exchange. For a lipoid insoluble small sugar molecule this time is of 
the order of 1-2 seconds. This is the actual average time spent in the 
capillaries. Thus we may use the anatomical figure for capillary length, 
call this length Z, of .7 mm. 

The volume supplied by a capillary is limited by the metabolic rate of 
the tissue. Consider this volume as consisting of a cylinder coaxial with the 
capillary. Consider the oxygen concentration at a distance x from the axis 
of the cylinder. This concentration is given in a formula presented by A. 
Krogh (1929, p. 270), 

To— = SPL 1.15? og 2 F"), (5) 
in which 7» and 7, are the oxygen pressures (in atmospheres) in the capil- 
lary and at the point x, respectively, d is the diffusion constant (number of 
cm.’ which in one minute will diffuse across 1 cm.” with a pressure gradient 
of 1 atm./u), p is the metabolic rate of tissue (number of cm.* of oxygen 
consumed/minute/cm.’ of tissue), R is the radius of the cylinder, and r is 
the radius of the capillary. 

To obtain the volume supplied by one capillary we desire the maximum 
distance 2) to which oxygen can be supplied under conditions of high levels 
of tissue oxygen consumption. This distance » will be evaluated later in 
this paper. 

III. The connecting system between the aorta and the capillaries. 

A. The structure of the branch system. This is the main problem of this 
paper. It is obvious that we desire some systematic way of splitting the 
single aorta into branches and do this in space so that the smallest 
branches (the capillaries) will regularly cover a given volume of space. 

An obvious method of splitting the large vessel is to split it into two 
branches; split each of these into two, etc. This is illustrated for three such 
splits in Figure 1. In general after 7 such splits there will be 2* small 
branches. 
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The problem of uniformly distributing these small branches over a 
given volume of space is more difficult. Rather than deal with the complex 
mammalian form as our volume of space, let us simplify the problem and 
deal with a simple abstract volume as a representation of the complex 
forms. Two such simple volumes are the cube and the sphere. Since the 
entire volume of the cube can be split into smaller and smaller cubes, 
which cannot be done with the sphere, we choose the cube as our abstract 
volume. Thus we have idealized our mammal into a cube. 

The obvious method of subdividing the cube into smaller cubes is by 
splitting it in eights by three planes parallel to the sides. This is illustrated 
for two successive splits in Figure 2. 


Kh A) 


FIGuRE 1 


Thus each split of the cube divides it into eight smaller units. However 
we want each split of the main vessel to give rise to two daughter tranekes 
This discrepancy can be resolved by considering each split of the cube as 
accompanied by three branchings of the vessel supplying the cube. For the 
sake of rigor we may have a systematic method of splitting each cube. 
This is illustrated in Figure 3. Figure 3A illustrates the first split by one 
plane parallel to one of the sides. The vessel courses to the center of the 
cube and there splits into two branches each one supplying one of the two 
rectangular volumes into which the original cube has been split. These 
two vessels then run to the centers of these volumes and there another 
split occurs. This is illustrated in Figure 3B. Figure 3C illustrates the 
process for the last of the three splits. The three together comprise the 
first division of the large cube into eight daughter cubes. Thus we have de- 
vised a method of systematically supplying the center of each of the eight 
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daughter cubes with a vessel. It is obvious that this method of division 
may be extended to successive splittings of the daughter cubes, etc. 

If we call the length of a side of the large cube S, after / splits the small- 
est cubes will be of (1/2')S, and of volume [(1/2").S]®. Since the large cube 
represents the metabolizing mammalian volume, we desire a sufficient 
number of splits of this cube so that a small cube volume is finally attained 
corresponding to the volume supplied by one capillary. The volume sup- 


original cube cube after ist split 


2nd split illustrated for one cube 
resulting from 1st split 


FIGURE 2 


plied by one capillary is +R?L, where L is the length of the cylinder sup- 
plied by one capillary and R is its radius. Thus we find that 


1 3 
[a s| See (6) 

B. Relative diameters of vessels at a branching. Consider the branching 
system of Figure 4. Let points A, B, C, and D be fixed with the distances 
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BC and BD equal. We want to determine the relative diameters of the 
branches AB and BC and BD which will minimize the total resistance of 
the system, if we assume Poiseuille flow and a fixed amount of material to 
be available for the pipe wall. We call /, the length of AB, / the length of 
BC and BD, 7, the radius of AB, 7, the radius of BC and BD, 6; the thick- 
ness of AB, and 6 the thickness of BC and BD. If R, is defined to be the 
total resistance of the system to flow, from the Poiseuille law of flow we 
find that 


Wet! (7) 


FIGURE 3A FIGURE 3B 


Ficure 3C 


where » and @ are defined as before to be the viscosity and density of the 
fluid, respectively. 
If M is the amount of material available for the walls of the entire 


system, then 
Mis 2a rhb t 4a rolods . (8) 


We must now determine in what manner the thickness of the vessel wall 
depends on the vessel radius. We know that the tension T in the wall obeys 
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the law 
L=pr, 
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(9) 


where # is the pressure in the vessel and r is the radius of the vessel. To 
keep uniform stress in the system it is convenient to have the vessel wall 
thickness depend directly on radius, a relation found to hold approxi- 


mately in the vessels of the animal. Thus 


6=ar, 
where a is a constant. 
From (10) we obtain 6; and & in terms of 7; and 7. Thus 


6;=ary and ds=aro. 


FIGURE 4 


Substituting (11) in (8) we obtain 
M=27ra ril, + 47a rl, 4 
or 
_ (M—2r0 ae 
Ae 4raly ; 


The condition of minimum resistance may be expressed as 


dR r 


wee 


But 
dRyr_ ORr OR df 


dr On Ors dr 


Substituting in (15) derivatives from (7) and (13) we obtain 


at 8n ols -2(—4 a) us ). 
mo 7 wo 7° De Ps 


(10) 


(11) 


(12) 


(13) 


(14) 


(15) 


(16) 
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This reduces to 


re —4 11 
ae aa 
or 
r= Ars. 
Thus we finally obtain 
f,=.1947;. (17) 


C. Relative lengths of the branches. Each split of the cube into eight 
smaller cubes corresponds to three branchings of the vessel supplying the 
large cube. Thus, since each split of the cube halves the linear dimensions 
of the component cubes, let every three branchings of the incoming vessel 
halve the lengths of the vessels. This can be accomplished in the following 


TABLE I 
Generation miei! Radius Leathe Resistance of single vessel 
of vessel thi aa ae of this generation 
ae is gener- =f; =]; =Ri 
i ation =n; try 
0 1 ro h (8/20) (lo/75) 
1 2 ro Io/21/8 (8n/mo)(Io/2/*)(1/y*70) 
2 + ro Iy/ 22/8 (8n/ mo) (bo/27/8)(1/-yr6) 
i 2s ¥'ro }y/2*/8 (8/20) (Ip /2*/8)(1/-y**r9) 


regular manner. If the initial vessel length is /o, then after i branchings the 
length would be J)/2*/, 

Before considering the total resistance of the system to flow let us as- 
semble the data we have so far obtained. This is done in Table I. 

IV. Total resistance of the system to flow. We consider the system as a 
series of many parallel vessels. After the ith branching we find 1; vessels, 
each with a resistance of Ri. Thus the series resistance of the vessels com- 
ing after the ith and before the i + 1st branching is (1/n,)Ri. 

The total resistance of the system, call it Rz, is thus 


> 1 1 1 
Rr=RotgRitERs...2-Ri--.5-Rm or Rr = DY——R:, (18) 
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where m is the total number of branchings considered. Substituting values 
for Ri from Table I in (18) we find 


ee, ee LE NAS 


Rr= 


8 
= n lS 1 (19) 


4 dad 241/38 41 
merge Meat) Seok 


Let 
1 
B 2S e2 (20) 
then 
¥ 8n lb fi— pm 
ei lies oe 


In (21) we have used the standard summation formula for the geometric 
series of (19) 

We are now in a position to calculate the resistance of the entire system 
once we know the parameters Jo, ro, y, and m. 

V. Comparison of our constructed system with actual anatomical system. 
Excellent data for the vascular system of the dog have been assembled by 
H. G. Green (1950). Because of the convenience of using his assembled 
data we have done our calculations for a 13 kg. dog. 

In the following paragraphs we will consider probable values for the 
parameters of the equations we have derived. 

The dog can be idealized into a cube of side of 23 cm. However, owing to 
the compactness of the model and the not necessarily direct paths of blood 
vessels it is convenient to choose Jp as 40 cm. 

The parameter 7) is the same as the aortic radius on which we have 
placed a limit in formula (4). In this equation C is the cardiac output 
which is given by Green as 2400 ml./min. (Joc. cit., p. 230) or 40 ml./sec., 
and n/a is the dynamic viscosity which is given as ranging from .017—.027 
cm.?/sec. (Joc. cit., p. 241), the high viscosity being used for the aorta and 
the lower for the small vessels. Substituting these values in (4) we find 

40 cm.3/sec. 


74 = (1100) x (.027 cm.*/sec.) (23) 


SAS Cie 


We have evaluated y on minimum resistance principles and found a 
value of y = .794 [see eq. (17)]. 
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The parameter m represents the total number of splittings of the vessel 
system. It may be evaluated most satisfactorily by considering the volume 
of tissue supplied by one capillary. Then consider the total number of such 
volumes comprising the entire animal, and from this we may easily obtain 
the number of splittings of one vessel necessary to achieve a number of 
vessels sufficient to supply this number of elementary volumes. The vol- 
ume supplied by one capillary is 7R?L, where, as in (6), L is the length of 
the cylinder supplied by one capillary and R is its radius. We have chosen 
Las.7 mm. and R may be evaluated from (5). Since we want point x to be 
at the outer surface of the cylinder, we let x = R. Thus we obtain 

104p 


y fa oP 
ea PU 1.15? log —* 7]. (23) 


Since the animal must be capable of sustained effort at high levels of 
metabolism, the largest volume supplied by one capillary must be evalu- 
ated from high oxygen consumption considerations. 

The average capillary concentration of oxygen is J and has a value of 
about 70 mm. Hg pressure or about .1 atm. pressure. The tension of 
oxygen at the surface of the cylinder is Tz, and since we are considering 
the maximum region that may be supplied, T, is 0 atm. 

The metabolic rate of the tissue (cm.*0:/min. cm.’ tissue) is p. This may 
be evaluated by considering the total oxygen consumption of the dog 
during rest (155 ml./min. from Green, 1950, p. 244), multiplying by a 
factor of 12 for the active state and dividing by the mass of the dog. Thus 

— (155 ml. 02/min.) (12) 
ms 13000 gm. 


p =.143 cm.? 0./ min./cm.} tissue . 


The diffusion constant is d and is given by R. Héber (1945, p. 273) as 
.12 cm.*0;/cm.? with a gradient of 1 atm./y. 

The capillary radius is r and is equal to 5u. 

Substituting the above parameter values in (23) we obtain 


104(.143) RRL Spy 
oles oe OR Te 
aa [ 1.15R log = i |. 


Solving for R we finally obtain 
R=35p. (24) 


With this value we find that the maximum cylindrical volume that can be 
supplied by one capillary is 


TRL = w (.1225 X 10-4 cm.?) (.07 cm.) } 


= 2.7X 10-8 cm3 oo 
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Substituting this volume in (6) and recognizing that each split J of the cube 
corresponds to three branchings of the vessel system, and substituting for 
S°, which is the volume of the cube 13,000 cm.%, we obtain 


[sea] S3 = 2.7X10-* em3 


or 
§3 


am — 
Sa Oa" Cin. s 


Then we find 
m=33. (26) 
We can now evaluate the total resistance of the system by substituting 
in (21), with 8 = 1 and n/o assuming an average value for the circulatory 
system of .022. Thus we obtain 


(022) cagya [34 | (27) 
See eC. SCC, Tk 


In obtaining (27) note that with 8 = 1 and summing fromi = 0 toi = 33 
in (19) there are 34 terms in the series and so its sum is 34. 
Total flow through the system can be obtained from the formula 
AP 
ie (28) 
where AP is the pressure head applied to the system, F is the total flow, 
and R, is defined as before. Substituting in (28) we find 


(10.0) (13) (980) ) 
2220 j (29) 


F = 57.5 gm./sec. 


ee 


In this formula we have used the physiological value of arterial pressure of 
the dog of 100 mm. Hg. The value of F obtained here is close to that 


observed in the dog (see p. 71). 
We may also calculate the radius of the capillary from our derived data 


of Table I: 


Teapillary = 133 

= 73375 
(.749) 33 (.43) cm. 
== 2. 2h 


(30) 


Green lists nine major types of vessels from the aorta to the capillaries. 
The thirty-three branches of our vessel system correspond to eleven major 
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splits of the original cube of the idealized dog. Thus there might seem to 
be some analogy between each split of the cube (or every three branchings 
of the vessels) and each new major type of vessel of Green’s classification. 
However, it is very difficult to make an accurate analogy of this type. 

We have accomplished a major part of our original aims. Utilizing a 
minimum number of parameters we have been able to arrive at satisfac- 
tory limiting values for important component parts of the vascular system. 
The total system we have constructed has given good agreement with 
anatomical data concerning vessel radii, lengths, and, most important of 
all, total resistance of the system giving good flow values. 

This investigation was supported by a research grant H-627 from The 
National Heart Institute, of the National Institutes of Health, Public 
Health Service. 
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Experimental procedures are suggested to test the theory developed in previous papers 
which may have applications to predicting the spread of information by contact through a 
population. The experiment is designed to test the statistical properties of the “acquaintance 
net’’ of the population. Thus behavioral aspects are for the time being eliminated, and atten- 
tion is focused exclusively in the properties of the potential communication net itself. 

Since the theory is not mathematically rigorous, it might be advisable to perform experi- 
ments with statistical materials only, such as playing cards, to test the validity of the assump- 
tions and approximations of the theory in dealing with particular statistical processes. One such 
experiment has been performed. The results agree closely with those calculated from one of 
the derived equations. 


In previous papers (Rapoport, 1953a; 1953b), hereafter called I and IT, 
we derived a number of models describing the spread of information by 
contact in a population, where instead of a thorough mixing a certain 
amount of ‘‘cohesion” occurred among the knowers. This bias is based, of 
course, on the well-known fact that the likely contacts of two individuals 
who are closely acquainted tend to be more overlapping than those of two 
arbitrarily selected individuals. 

The “assumption of transitivity’’ in its strongest form (cf. I) led to a 
spread of state drastically slower than that derived from the assump- 
tion of randomness. In fact, the spread on the basis of that assumption 
was considerably slower than that observed experimentally, which, in 
turn, was considerably slower than the spread in a random net. 

These observations led us to introduce modifying parameters in various 
ways to “weaken” the assumption of transitivity, so as to get situations 
intermediate between completely “‘transitive” and completely random 
cases. In particular, two forms of the ordinal time course equation gave 
excellent agreement with experimental results, namely, 


e(t+1) =1- [1-2] [1-P() (1— er) J emo (1) 
75 
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and 


e(¢+1)=f— [f—2@)) (2) 


As in previous papers (é) is the total fraction of knowers, and P(#) the 
fraction of new knowers on the ¢th remove. Wherever the argument is 
omitted, it is understood to be t. 

The parameter a in both equations is the “actual axon density,” Le., 
the total average number of contacts per individual during the entire 
process. In equation (1) the parameter 6 measures the inclination to seek 
contacts outside one’s acquaintance circle as the saturation of the popula- 
tion with knowers increases (i.e., 8 measures the motivation to spread the 
information as far as possible). In equation (2), the partial randomness of 
the process is already “built in” by the assumption of equidistribution of 
new knowers among the group of knowers. The parameter fin that equation 
measures a different quantity, namely, the fraction of the population who 
participate in the process at all. 

In spite of the different assumptions underlying equations (1) and (2), 
they lead to very closely coincident ordinal time courses of the spread so 
that in the given situation it is impossible to choose between them. The 
equations do indicate, however, how an experiment might be designed to 
test either hypothesis. Our whole theory of the spread of information by 
“removes” is based on the theoretical tracing of the communication net 
of the population. It is assumed that information will-flow along these 
channels. It is, however, possible to test our theory without this assump- 
tion, namely, to test the statistical properties of the communication net 
itself. We specifically equate the communication net with the acquaint- 
ance net, well realizing that for different situations the net may be dif- 
ferent (for example, an anecdote may travel by channels entirely different 
from those traveled by a rumor, and various rumors may travel by di- 
verse channels). 

At any rate, we propose to apply our theories to an experimentally de- 
termined acquaintance net. For this purpose, it is necessary to take a 
single poll of a population. Let each individual of a large population, say, 
a university faculty of some 10® members, name 10 other members in the 
approximate order of frequency of his contacts with them. Let us assume 
that we have this information from all the members. We are then in pos- 
session of a thousand cards, each card having on it ten names. The rest 
of the experiment can be performed without any further contact with the 
subjects. 

We select a certain fraction of the cards at random and call that frac- 
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tion P(0). Next we fix a arbitrarily. For example, let a = 2. We then note 
the first two names on each of our initially selected cards and pull all the 
corresponding cards out. Those not included among the P(0) are the in- 
dividuals of P(1). We repeat the process with the first two names on the 
cards of P(1), etc., thus obtaining the whole set of P’s and therefore of 
«x’s. We then proceed to test the applicability of equation (1). Since the 
“motivation” reflected in the parameter @ is absent in our experiment, it 
should not appear in the equation. However, the quantity 6, which in 
equation (1) appears as e~**, need not be unity. There may be a constant 
amount of randomness in that process. Therefore, let us keep the parame- 
ter 6 as a constant <1. Equation (1) then becomes 


Perea ee ey (2) (1 et) | eer (3) 


for ¢ > 1. If the first set of P’s are well fitted by equation (3) for some 
value of 6 (0 < @ < 1), we can put the underlying theory to much more 
severe tests in many different ways. 

1. Repeat the experiment an arbitrary number of times for different 
sets of initially selected cards, keeping the number of these cards constant. 
If the theory is correct, the same value of 6 (which is supposed to be 
characteristic of the population) should fit all the sets of data. This is a 
test of the significance of @ as an invariant characteristic of the popula- 
tion. 

2. Repeat the experiment for different values of P(0). Here we deliber- 
ately vary one parameter to see if the previously established values of the 
other parameters can still fit the new sets of data. 

3. Repeat the experiment for different values of a, getting still another 
family of curves, which, if the theory is correct, ought to be fitted by the 
same values of the other parameters. N.B.: the parameter a need not be 
confined to integral values. We can, for example, let a = } by choosing 
alternately 1 and 2 names on successive cards. 

Next we can test the significance of the parameter @. As we have seen, 
@ is a measure of randomness in the acquaintance structure, or rather the 
lack of it, since 6 = 0 for completely random structures. Keeping a con- 
stant, therefore, we select instead of the first a names on each card, a 
names further down the list. If the position on the list is roughly an indica- 
tion of the “tightness” of the acquaintance relation, it should follow that 
the further down on the list the names one selects, the smaller should be 
the @ which fits the corresponding set of data. Thus the significance of 0 
could be tested at least qualitatively. 
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The choice between equations (3) and (2) to fit the set of data can be 
decided by the behavior of the apparent axon density 
ares 1— x(t) (4) 
a) =p iy Toe G+ I)’ 
We have seen in II that if equation (3) determines the process, we should 
have 
pis 
P (t) 
for t > 1. That is to say, a(é) varies weakly with P(#), and for small P() 
remains approximately constant and equal to 


1—e-%+ (1-—4)a, for yeas Ale (6) 


a(t) = log Li = Pt) merOh athe (S) 


Therefore the apparent constancy of a(¢) for > 1 should indicate the use 
of equation (3). However a(1) may be different from a(0), and the size 
of 6 is an indication of this difference. 
On the other hand, equation (2) implies that* 
pet =f — e—4P/z >1 7 
alt) —57y logit at i= & ve fe C8) 
Expanding the exponential and logarithm in (7) we can write, as a first 
approximation, 
a*P 
a(t) Yya-—-—, (8) 
2% 


which becomes increasingly better as P grows smaller and x grows larger. 
But x, being cumulative, is monotone increasing from the start, and P is 
eventually monotone decreasing. Therefore (7) yields an eventually in- 
creasing apparent axon density. If this is observed in the data, this is an 
indication to use equation (2). 

In equation (2), the parameter f denotes the fraction of the population 
participating. This parameter can be introduced into our experiment by 
arbitrarily marking a fraction of the cards “non-cooperators.” The param- 
eter f can be arbitrarily varied and so a family of experimental curves 
obtained. 

One other possibility should be mentioned, namely, that @ in equation 
(3) may depend on a. This is not apparent in the formulation of the theory, 
but is apparent from a qualitative argument. If a is small (say, between 
1 and 2), the probability that the process is random should be greater 


* If f represents the fraction of participants, the generalized definition of a(t) becomes 


5 Shoe AU, 


— 


PO) 8 F—e0Fiy- 
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than if @ is so large as to include a majority of one’s acquaintances. There- 
fore we may take a function of a which is zero for vanishing a and 1 (or 
some other constant) for large a. We might take, for example, 


Cie (9) 


This leads to another possible form of equation (3), namely, 
a(t-+1) =1— [1—-«] [1—P(#) (1— e-2-e*)) | e-2e*P | (10) 


This form is significant if we perform the experiment for several distinct 
values of a. 

It should be noted that the proposed experiment is not designed to test 
the actual dynamics of the spread of information through a population 
but only the statistical aspects of an acquaintance net in terms of our 
random net theory modified by the ‘“‘assumption of transitivity” and cer- 
tain other parameters. In a way, such an experiment would be a pre- 
liminary investigation. If the statistical behavior of the acquaintance net 
is shown to be adequately described by one of our equations, it might be 
worthwhile to design further experiments to see whether the dynamics of 
information spread follows similar laws. 

It seems advisable to perform experiments such as the one suggested 
for still another reason. Our probabilistic reasoning, in I and II, as will 
appear to any critical reader, was far from rigorous. Indeed rather arbi- 
trary assumptions were introduced, simply because a completely rigorous 
formulation was beyond the present capacities of the author. That is to 
say, one might perform experiments simply to test the applicability of the 
equations in some probabilistic process. Numerous physical models have 
been designed to trace “artificial epidemics.” As far as the writer knows, 
a thorough mixing of the population was always assumed in these models. 
Our theory, however, is specifically designed to describe a process without 
thorough mixing. The following model was constructed to test just such 
a theory. 

Spread of a state through a deck of playing cards without shuffling. We 
take two identical decks of cards, one of which will serve as our “popula- 
tion,” and the other as an image of the population to be used only for 
convenience of keeping records. We arrange the first deck as a single row 
or a circle, having first selected two cards at random as the starters, which 
are then removed from the second deck and called P(0). The experiment 
can be performed quite quickly if the second deck is arranged in order, so 
that the new knowers, indicated by observing the first deck, can be quickly 
found in it. 
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Now we suppose that the four neighbors of each card, two on the left 
and two on the right, are the only acquaintances, so that only to these can 
information be passed. These four are labeled 1-4 from left to right. At 
each remove, each new knower passes the information to randomly se- 
lected two of his four acquaintances. Which two these shall be can be de- 
cided by some random procedure. At each remove, a new group of new 
knowers is taken from the second deck. Then lots are cast to decide for 
each of them to which two of his four acquaintances he will tell, and these 
in turn are removed from the second deck if they had not already been 
removed. 

In the experiment conducted by the writer, a sequence of numbers con- 
sisting of 1’s, 2’s, 3’s, and 4’s was determined by the appearance of the 
four suits in drawing from a well-shuffled deck. Besides these another such 


TABLE I 


t 0 1 2 3 4 5 6 | 7 8 | 9 


“(t)obs:| 2038 |. 103k 177 e258" | a2aO ole 1 aooOnt S400 coos Me soo 


x(t) calc} .038 | .103 | .166 | .222 | .269 | .306 | .334 | .345 | .368 | .376 


sequence was written, and the pairs so determined constituted the de- 
cisions to which acquaintance each knower must tell. Of course a pair 
(i, 7) meant that only one acquaintance would be told. Since the prob- 
ability of (2, 7) is +, this would mean that in an infinite population the 
parameter a would have the value 1.75. In our population of 53 cards (in- 
cluding joker), we would expect a to be somewhat greater. In fact, as de- 
termined from the data, a turned out to be 1.80. 

If the cards are not shuffled, a strong bias would be expected to work, 
and possibly the “assumption of transitivity” could be applied. As it 
turned out, the data could be fitted very well for 6 = 0.8. The average 
results of ten runs are shown in Table I. As is apparent, a(é) after the 
usual initial drop remains approximately constant. Therefore the use of 
equation (3) is indicated. By letting 6 = 0.82, we obtain a(é) = 1.1 for 
¢ > 1 and the calculated values of «(#) as shown in Table I. 

It is to be expected that if shuffling is introduced, either the value of 
8 will become lower or 4 will not be constant but a decreasing function of 
t (decreasing with each shuffling). By introducing some standard shuffling 
procedures, the ‘“‘amount”’ of shuffling can be measured in some way and 
the relation between this amount and the value of @ or one of the other 
analogous parameters can be determined. 
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A systematic experimental procedure of this sort should establish one 
or the other of our equations (possibly with some further modifications) 
and furthermore give more precise meaning to the parameters. Once this 
“phenomenological” theory is sufficiently developed to be of some value 
in predicting the outcomes of such purely statistical experiments, it might 
be advisable to carry the theory into the field of social statistics to see 
where it might apply. 

This investigation is part of the work done under Contract No. AF 
19(122)-161 between the U.S. Air Force Cambridge Research Labora- 
tories and the University of Chicago. 
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Onsager’s method of studying the nematic phase is developed for general molecular interac- 
tions. It is shown that the symmetry of the molecule helps determine the type of transition 
that occurs in passing from the isotropic phase to the anisotropic phase. 

The possible relation between the nematic phase and spindle structure is briefly discussed. 


I. Introduction. Studies of the structure of spindle and the changes 
that occur in the spindle during the course of mitosis and meiosis provide 
a powerful stimulus to theoretical and experimental investigations of the 
mesomorphic, or “liquid-crystal,” states of aggregation. The various 
parts of the spindle have a number of properties in common with the 
mesomorphic states; structural rigidity and flow, optical birefringence, 
changes of consistency and birefringence with changes in concentration 
of either macromolecules or smaller molecules all lead one to conclude 
that the various parts of the spindle are either liquid crystals themselves 
or else are so remarkably analogous that a study of liquid crystals may 
throw some light on the structure of the spindle. Thus J. D. Bernal 
(1940) has suggested that the spindle may be tactoidal in nature. G. Oster- 
gren (1949) has emphasized the “fluid” nature of spindle fibers. 

While a well known fact, it should be emphasized that the spindle is 
never homogeneous. This is clearly demonstrated, for example, in the work 
on birefringence of M. M. Swann (1952) and the work of J. G. Carlson 
(1952) on microdissection. Furthermore the properties of various parts of 
the spindle change with time, being closely correlated with the position 
and motion of the chromosomes. Therefore one cannot deal with the 
spindle as a liquid crystal. One can only treat the various parts of the 
spindle as being in or analogous to a mesomorphic state. 

The tactoid has received the most attention in the literature (see, for 
example, Schrader, 1944) as being a state of matter closely resembling 

-what is known about spindles. However the tactoidal state is only one of 
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many states having the requisite properties and a systematic attempt at 
investigating the other mesomorphic states for possible applications to 
spindle structure should be attempted. One of the difficulties in such an 
attempt is the scarcity of statistical mechanical theories of the various 
types of liquid crystals. 

IT. Classification of mesomor phic states. G. Friedel and E. Friedel (1931) 
classify mesomorphic states as nematic or smectic and, of the nematic 
type, they recognize two kinds—the true nematic and the cholesteric. 
The nematic type is distinguished by the presence of long chain molecules 
or long rod-like colloidal particles which are more or less aligned in a par- 
ticular direction in some solution or matrix. Furthermore the centers of 
gravity of the molecules are randomly arranged, i.e., they show no 
periodicity in any direction. The cholesteric differs from the true nematic 
in that the cholesteric shows strong optical activity while the true nematic 
does not. Presumably, in the cholesteric there exists a screw-like arrange- 
ment of important elements of the molecule. 

The smectic state also contains rod-like elements which are more or 
less aligned. Here, however, the centers of gravity of the elements are not 
randomly located. Along the direction of orientation the centers of gravity 
form a one-dimensional periodic lattice; however, in the plane perpendicu- 
lar to this direction the molecules are oriented at random. 

To the types just described there must be added the tactoidal structure. 
Here there is no periodicity along the direction of orientation, while in 
the plane perpendicular to this direction the molecules or colloidal par- 
ticles form a hexagonal lattice. 

In brief, in the nematic state there is no periodicity; in the smectic 
state there is one-dimensional periodicity, while in the tactoidal state 
there is two-dimensional periodicity. All the mesomorphic states have the 
common property of demonstrating alignment along some direction in 
space; thus all of them exhibit birefringence. 

Substances in the mesomorphic state may vary considerably in con- 
sistency, ranging anywhere from a thick paste to a thin watery solution. 
Furthermore, the distance between molecules, and, hence, the consist- 
ency, is quite sensitive to the ionic concentration of the medium and the 
temperature; even the existence of any particular mesomorphic state ne- 
cessitates having the specimen between certain limits of ionic concentra- 
tion and also within a well defined temperature range. 

All substances in a mesomorphic state exhibit elastic properties under 
low stress, but yield at higher stress and show flow properties. This phe- 
nomenon alone would make them of interest to cytologists and cell physi- 
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ologists as examples of solutions having deformation and flow properties 
similar to protoplasm. The biologic interest, however, is heightened by 
the many properties which such substances have in common with por- 
tions of the mitotic spindle. 

III, The nematic phase. No claim can be made that any part of the 
spindle is actually in a nematic state. The submicroscopic structure of the 
spindle is as yet unknown. However the nematic state appears to be the 
simplest from a theoretical standpoint, and it appears best to investigate 
this state further before studying the others. Furthermore, it may turn 
out that certain general properties of oriented structures will suggest 
themselves as a result of a study of this special case. 

For our purposes the nematic phase will be defined as a solution of 
particles containing an axis of symmetry, so arranged that the particles 
exhibit no periodicity. 

The method used here will be that of L. Onsager (1949) and a part of 
his analysis will be duplicated in the interest of continuity and clarity. 

Let us now consider a system containing particles with an axis of sym- 
metry, say, for definiteness, long chain molecules. The system will also 
contain other particles, but our attention will be focussed on the former 
particles only. We shall adopt Onsager’s terminology, which is one that 
reflects our interest in these particles, i.e., we shall call them ¢he particles 
and shall refer to the rest of the system as the solvent. The entire system, 
i.e., the solvent plus the particles, will be designated as the solution. 

Let w be the potential of the forces between the particles, averaged over 
all positions of the other molecules of the solution (Fowler and Guggen- 
heim, 1949, p. 259 ff.). 

Let F = excess free energy of the solution over that of the solvent. 

N = number of particles in the volume V. 
Z = configuration integral 


== pain femtral sg? ocd ys 


P = osmotic pressure of the solution. 
yw = chemical potential of the particles. 
Then we have the thermodynamic relations 
F= Nwo— kT logZ, (1) 
where p, is a constant depending only on the temperature. 
Also 
an é | 
fps ir| 7 (log Z) 


T:N 


86 IRVIN ISENBERG 


and 
oer ( (log Z) : (2) 


TsV 


Onsager introduces an angular distribution function f(a) which may 
be defined by stating that f(@)dQ gives the fraction of all the particles 
which point into the element of solid angle dQ surrounding da. Thus the 
number of particles dV pointing into dQ about a is given by 


dN=Nf(a)dQ. (3) 


It is clear that f must be normalized, i.e., 


Sf(@)da=1. (4) 
Furthermore it is clear that if the solution is isotropic, then 
1 
= =. 5 
f = constant Ag CS) 


Suppose now that one considers an arbitrary f(@), i.e., one which is 
not necessarily that which an equilibrium solution at some fixed T and V 
would have. Then the free energy corresponding to this arbitrary f(@) 
would be higher than the equilibrium free energy, or, what is the same 
thing, log Z would be lower. The Onsager method consists, in part, of a 
calculation of log Z for an arbitrary f(a) under the condition that the 
solution is dilute, a condition which is realized in many substances in a 
mesomorphic state. He also assumes a particular form for w. The reader 
is referred to Onsager’s paper for a discussion of the form of w which is 
chosen. In this paper w will be left quite general. This generality has sev- 
eral advantages, two of which stand out. First, it should be noted that the 
precise nature of the interaction between the particles is not known. 
Second, as will be shown later, leaving w general permits the elucidation 
of an interesting relation between the symmetry of the particles and the 
nature of the isotropic-nematic phase shift. 

Under the conditions described above, Onsager shows that for arbi- 
trary f(a) 


loe Z = yn} 1 —log PSs @) log [4a f (a) ] dQ 


+558 (4, a’) f(a) ¢ @’) dada’ +...|, 
where 


B=TSdVifd Va lew? — 1] (7) 
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and w is here the interaction between two particles whose orientation is 
fixed. The integrals in equation (7) are over all possible positions of the 
centers of gravity of the particles, the orientations remaining constant. 
Let us assume that w depends only on 6, the angle between the particles, 
and the distance between the centers of gravity. Then 8 will be a function 
of @ with 


b=cos 6=a-a’. (8) 


As indicated in (6), there are neglected terms on the right-hand side 
of this equation. These are terms with coefficients (V/V)?, (V/V), etc. 
Under the assumption that the solution is dilute, these terms will be neg- 
lected. 


The following notation will be introduced 


c = concentration of particles = + (9) 
| f| =Sf @) log [4x f (2) ] da (10) 
gl fl = —SSB(u) f(a) F(a’) dada’ (11) 
p=ots cg. (12) 


The notation illustrated in equations (10) and (11) is designed to show 
explicitly that o and g are functionals of f. 

The problem now facing us is to find the equilibrium distribution func- 
tion f(@) by minimizing y subject to the restriction 


ffdQ=1. 
This is a variation problem. If we let v be a Lagrangian multiplier we have 
do|f| +2 cég|f| —vdffda=0, (13) 
which leads to 
log 4a f +1—cfB(u) f(a’) do’ —v=0. (14) 
Equation (14) is an integral equation for f. As will become evident, the 


solution of (14) is not unique. For example, it is always solved by 


1 
f = constant = ae 
since 


SB (pw) da’ 


is independent of a. 
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For the isotropic distribution 


1 
7|7-| =9. 


As shown in the Appendix, o is positive for any other distribution. It is 
also shown there that o has many properties which make it a natural 
measure of the degree of orientation of the system. A system more highly 
oriented will have a larger o value than one less highly oriented. Therefore, 


it may be considered appropriate to call o the order of the system. 


Consider the following variation problem which will be referred to in 
the remainder of this paper as the modified variation problem. Here the 


value of o|f| is fixed, say, 


g | f\ mee 
o being a given number. We still demand that 
f7fdQ=1, 


and we wish to minimize g|/|. 
This problem leads to the integral equation 


2 2v 
bg+5 b0—> sf fda=0, 


where 2/k and 2/k are Lagrangian multipliers. 
It will be seen that (15) is the same as (14). 
Where a solution of equation (15) reads 


f= f(a; ky»), 
the corresponding solution will be 
f= f(G; ¢,y). 


(15) 


(16) 


(17) 


The new variation problem attempts to solve the previous minimum 
problem in two steps. Instead of minimizing y directly, ¢ is held constant 
while y is minimized by minimizing g. The next step is to minimize y with 


respect to the value of ¢ originally chosen. 
The condition 
ffdQ=1 


establishes a relation between & and v, and we may therefore consider 


f= f(a; B). 


Furthermore, the definition of g and o and the restriction 


o|fl=o, 
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where o is a given value, yields equations 
Gr—t04 0K) (18) 


and 
g=g(k). (19) 


These last two equations may be considered parametric equations for g as 
a function of o. 

What we have done is to fix the order and choose that distribution 
which makes g a minimum. What we must now do is to fix the concentra- 
tion and choose that order which makes y a minimum. This is equivalent 
to choosing that value of & which makes 


¥(k) =o(k) +3 cg(k) (20) 
a minimum. 
We then have 
do (k) dg(k) _ 
<7, t+} 0-4 =0 (21) 
d?a(k) ,, _ d?g(h) 


hy ae a Sar pe >0. (22) 


Equations which are equivalent to the two above are 


# tots cg(o)1 =0, 


d2 
a Aes cg(c)] >0, 
or 4 
ite 0 (23) 
o 
and - 
aS 24 
= pee (24) 


The distribution function f(@) may be expanded in a series of Legendre 
polynomials 


4a f(a) = Si (21+1) A:Pi(cos 6) , (25) 
1=0 


where A, = 1 to provide for the proper normalization of f(a). 
The function B() may likewise be so expanded 


Bi) = > DP) (26) 
. +=0 . 
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It is then not difficult to show by the use of the addition theorem for 
Legendre functions (Copson, 1935, p. 304) that 


glfl=— > Dsl. (27) 


Likewise it may also be demonstrated that 


e\sl=520 (27+-1) Ai—28 42—-Fida— --- +GiAi sgh 
I=1 


+6,At+G,Ait+...—3AjA2—9A1A2As—... - 


The numerical constants F; and G; are all positive. Their precise value 
is of no interest to us. Several things should be noted about (28). The 
cube terms such as A3, A, etc., all have negative coefficients and involve 
only terms with even subscripts. The terms of the fourth power are all 
positive and all subscripts, both odd and even, are represented. 

The modified variation problem may now be attempted, using the 
representations of equations (27) and (28). The restriction that f(@) be 
normalized has been previously taken into account by equation (25). 
Since we are now dealing only with distribution functions that are nor- 
malized, the modified variation problem will read 


sel sl+< sol s/=0, (29) 


where 2/) is a Lagrangian multiplier and the variation of the distribution 
function is to be accomplished by varying all the 4,’s. 
Equations (27), (28), and (29) lead to the set of equations 


(3 — Dj) A,;—6A,;A.—9A2A3+. a 8 =a) 
(S —AD2) d2a—94143—28 Apt... =0. 


(30) 


These equations may be viewed as being equations for the A,’s, a solu- 
tion of which will yield A; in terms of \. It can be seen that A; = 0 for 
all 7 constitutes a solution to the set. This is the isotropic solution that we 
have previously seen solves the variation problem. 

It must be emphasized that we are here dealing with the modified varia- 
tion problem; the value of o in equation (28) is preassigned. If the tem- 
perature, osmotic pressure, and concentration of a solution are given, 
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then we know, of course, that the assigned value of o will not, in general, 
correspond to the actual value of «. However, from the manner in which 
the Lagrangian multiplier was chosen we know that when z is the actual 
concentration ¢ will be the actual order of the system. 

It should also be realized that \, defined by (29), is the same as the 
parameter & of equations (18) and (19). 

For reasons which will become clearer later we shall be interested in the 
behavior of the set of equations (30) in the neighborhood of o = 0, or, 
what is equivalent, where all A,’s are very small. In this region we have 


(SN Ag 0 


(S —AD2) Ar=0 
(31) 


Suppose now that \ is so small that all the brackets [(27 + 1) — XD] 
are positive. Then the only solution is A; = 0 for all 7. Now imagine that 


Oo P » 


FIGurRE 1 


\ is made larger until one of the brackets vanishes, say, the one with in- 
dex n. At this point there are two solutions, the new solution being one 
for which 


A,;=0, 1An, \ (32) 


A,#0. 


This new solution is an anisotropic solution to the set (31). To visualize 
this more clearly let us imagine a plot of o vs. \ (Fig. 1). 

From 0 to P there is but one value of o, namely, o = 0. For \ > P, 
a is double valued; o may still be identically zero or o may rise with 2. 

We wish to examine o and g in the neighborhood of P, i.e., in the neigh- 
borhood of \ = (2% + 1)/(D,). In particular we shall concentrate on the 
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otropic solution, the properties of the isotropic solution being evident. 


anis 
Two cases must be distinguished depending on whether k is even or odd 
4n f (@) =1+ (2k+1) AxPx (cos 9) (33) 
g| f| = —DrAi— Dodo (34) 
2 
o|f|= PH" AL-PAL+ .. .if Bis even (35) 
2k+1 
o|f|= = Az+GiAi+. -.if kis odd (36) 
cs 
Cy 
es 
Coa 


FIGURE 2 


For later use the first and second derivatives of g with respect to oa, 
evaluated at A, = 0, will now be computed 


Bi 5 
Tope ge |e eee ue 
dA, J4,= 
d*g cavans 
_ 4) _.+ o if k is even (38) 
= d*g ws 16G;,D, . . 


It is shown in the Appendix that —dg/do is always positive. Further- 
more, from (37) it is clear that —dg/do must be finite at ¢ = 0. Thus, for 
example, —dg/do vs. ¢ might appear as shown in Figure 2. Also Shows in 
Figure 2 are plots of 2/c which are, of course, horizontal straight lines 
Equation (23) shows that an equilibrium point is found by noting the 
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intersection of the line 2/c with —dg/dc. It is clear that for sufficiently 
low concentrations there is no such intersection for the curve of Figure 2. 
When an intersection does occur, there are two points of intersection. 
However it is important to realize that the point P is unstable because of 
the stability requirement [eq. (24)] that d?g/do? > 0. It is therefore clear 
that only points to the right of R are stable and therefore that only aniso- 
tropic solutions with a finite degree of order will occur. As a consequence the 
transition from the isotropic solution to the anisotropic solution must be dis- 
continuous in the situation illustrated by Figure 2. 

The significant fact of the plot shown in Figure 2 is that it was not ar- 
bitrarily selected. It corresponds to the case of even k previously con- 
sidered. That this is so in the neighborhood of ¢ = 0 is shown by (37) 
and (38). The only remaining point that needs demonstration is that a 
maximum of —dg/do exists. But a maximum must exist, for if it did not 
—dg/do would rise indefinitely and the stability condition, (24), could 
never be achieved. A nematic state would never appear. Since such a 
situation is eliminated by the obvious a priori reason that we are develop- 
ing a theory of the nematic state and of the transition to the nematic 
state, it is clear that we must require —dg/do to have a maximum. 

It is, of course, conceivable that —dg/do should have more than one 
maximum. If such were the case, the possibility of a discontinuous transi- 
tion from one nematic state to another would arise. There appears to be 
nothing in the theoretical considerations presented here to rule out sucha 
possibility. The existence of more than one nematic state and, for that 
matter, more than one smectic state, has been questioned by G. and E. 
Friedel (loc. cit.). However, D. Vorlander (1933) has discovered substances 
with as many as four mesomorphic phases, and the possibility exists 
that more than one nematic phase may actually occur. The difficulty in 
visualizing more than one nematic phase may perhaps be due to the error 
of thinking of a nematic phase as containing particles all perfectly aligned. 
Actually there is, of course, no perfect alignment. Perfect alignment 
would correspond to ¢ = o, whereas all actual solutions have a finite o. 
A transition from one nematic phase to another would be a transition to a 
phase with a higher a, i.e., a transition to a phase with a higher degree of 
alignment. 

We now turn to the case of odd &. Here, recalling (39), it is clear that 
we are dealing with a plot such as shown in Figure 3. In such a case all 
points on —dg/do are stable and the transition to the nematic state 1s con- 
tinuous. Thus, depending on whether & is even or odd, the transition will 
be discontinuous or continuous. 
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It is possible of course that —dg/do in the odd case will have a maxi- 
mum somewhere that will be higher than 


In such a case the transition will be discontinuous. However, without 
eliminating this possibility, in the future the case of even & will be identi- 
fied with a discontinuous transition, while that of odd & will be treated as 
involving a continuous transition. 

Thus far hardly anything has been said concerning w, the interaction 
between two particles. The interaction, of course, determines everything 
about the nematic phase and it is important to see the relation of the 


ie 
Cy 


Q|N 


FIGURE 3 


interaction to what has been developed by way of a theory of the nematic 
phase. In particular we shall try, insofar as it appears possible, to under- 
stand the distinction between discontinuous and continuous transitions 
in terms of the interaction between molecules. 

From equations (7) and (26) we have 


1| bee} 
pla Visa Valew*r — 1] = Dale D:P;(u) . (40) 


If, in the left-hand side, we choose a coordinate axis with the z axis pass- 
ing through the line of symmetry of one particle, we may write 


Slew*? —1] dV = S" DP. (u) . (41) 


i=0 


The first thing we note is that if w(«) = w(—y), ie., if the interaction 
does not change when one particle is rotated through 180° about an axis 
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perpendicular to it, then D; = 0 if is odd, and, therefore, only even terms 
appear in the Legendre series. In this case the transition is discontinuous. 
This situation will prevail if the particles are long chain molecules with a 
center of symmetry. On the other hand, if there is a high degree of asym- 
metry then the largest D; will be odd and therefore the smallest quantity 
27 + 1/D; will probably correspond to an odd i. 

Therefore, roughly speaking, symmetrical interactions will lead to dis- 
continuous transitions while asymmetrical interactions will lead to con- 
tinuous transitions. 


APPENDIX 


The functional o|f| has the form of a negative entropy and, therefore, 
the demonstrations of properties II and III are identical to demonstra- 


tions given for the corresponding properties of entropy (see, for example, 
Tolman, 1938, p. 165 ff.), 


o|f|=JSflog(4af)dQ, where f>0 and ffdQ=1. 
Property I: 


Property IT: : 


Property III: Given any f = f, it is always possible to choose an f = fe 
such that fo fluctuates rapidly about f; and such that 


o| fol >o| fil 
glfel=elfil - 


We are now in a position to prove that 


$8 <9. 


o 
Suppose we have already found an f such that 
SfdQa=1, 


o| f| = 01 (given) , 
and - 
g|f| =minimum . 
Now let us choose a slightly higher value of oi. We can, using property 
III, pick a new f which will give the same value of g, but which will give 
the new and higher value of o1. But this choice of f has yielded a value of g 
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which is higher than (or equal to) the new {minimum- Consequently, by 
choosing a higher value of o; we get a lower (or equal) value of g 


wile 
tga 0: 
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Some conclusions of irreversible thermodynamics are summarized. It is shown that @, the 
rate of irreversible entropy production, is not minimized in the steady state. It is also postu- 
lated that multiple steady states are possible in nonlinear kinetic systems, giving rise to situa- 
tions of possible biological interest. The necessity of examining particular kinetic models is 
mentioned. 


For some time now attention has been drawn to the fact that all living 
systems are “‘open systems” and are usually in the “steady state.” As a 
consequence, various attempts have been made to deduce some of the 
properties of life from thermodynamic and kinetic considerations of open 
and steady state systems. Probably the best example of such an attempt 
is the work of A. C. Burton (1939). 

The purpose of this note is to summarize some of the conclusions 
reached concerning these properties and to correct what appear to this 
author to be some misconceptions in the literature. 

It seems worthwhile to begin with a number of definitions, since some 
of the misconceptions have arisen from ambiguous terminology. 

We shall call an isolated system a system exchanging neither heat nor 
matter with the environment. A closed system is a system exchanging heat 
but no matter, and an open system a system exchanging both heat and 
matter with the environment. Furthermore, we shall denote the flux of the 
ith substance across a unit surface normal to the flux vector as /;. In the 
case of a conserved incompressible substance the continuity equation 
states that 

div i+ 54=0, (1) 


where c; is the concentration of the ith substance inside the system. 
In the case of a non-conserved substance, i.e., if sources or sinks are 
present, (1) extends to 
div ct = 35, (2) 
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where q; represents the rate of production or consumption of the ith sub- 
stance. In the case in which J; is —grad c;, as, for example, in the flow of 
heat or diffusion, (2) becomes 


0¢; 


=V2c; ns 3 
Pe Vics 4 (3) 


A system is in equilibrium if, for all J;, JF; = 0. A system is in guasi- 
equilibrium if, for alli = 1, 2,...k, F; is a non-zero constant and, for 
alli=k+41,...n, f;: = 0. A system is in the steady state if, for all 
i, dc;/t = 0, or, according to (2), div Fi + qi = 0. 

We will summarize some of the conclusions of thermodynamics con- 
cerning open systems. It is well known (De Groot, 1951) that modern 
thermodynamics has changed the inequality of classical thermodynamics, 
dS > dq/T to an equality, dS = dq/T + dgq*/T, where dq* is the so-called 
“uncompensated heat.” It has turned out that this new term is more than 
just a formal definition. It has been possible to establish its value in many 
different irreversible processes (De Groot, Joc. cit.; Denbigh, 1951). One 
limitation of this method is that the end points of the process in question 
have to be equilibrium points (Tolman and Fine, 1948). But, more re- 
cently, P. W. Bridgman (1950) has proposed a “generalized entropy” 
whereby even this limitation can be abandoned. 

Another way in which to see the modern conception of entropy is to 
“concretize” entropy as if it were some kind of a non-conserved substance 
and then write down a continuity equation. We could then try to write 

a = 6+div S, (4) 
where @ is the rate of irreversible entropy production and div § the entropy 
flux across a unit surface. It turns out from more detailed physical con- 
siderations that equation (4) is actually justified. In evaluating its mean- 
ing, we must remember that the second law of thermodynamics requires 
that 6 be positive for all irreversible processes. Hence, given an open sys- 
tem (i.e., div § not zero), dS/dé can increase, decrease, or remain constant 
depending on the relative values of div $ and @. This is therefore the first 
generalization which can be made concerning open systems: The rate of 
change of their entropy is in no way restricted by the second law of 
thermodynamics. 

It has been stated by some authors (Prigogine and Wiame, 1946; De 
Groot, loc. cit.; von Bertalanfly, 1950) that the steady state is char- 
acterized by a minimization of @. Actually what has been shown is that 6 
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is minimized for the quasi-equilibrium state. This can be seen as follows. 
In the theory of irreversible processes, 70 can be written 


T= SFX, (5) 


where X; is the ith thermodynamic force (e.g., temperature gradient). If 

we now assume that all 7; can be written as linear sums of the forces, 

ji= > L,;X;, it follows by differentiation with respect to X,, keeping 
. 


all other X; constant, that 


BH 7 2 abs he. (6) 

Therefore 6 is minimized for all 7; = 0, a situation we have defined 
before as quasi-equilibrium. It is also worthwhile to point out with Den- 
bigh (loc. cit.) that (6) holds only if 7; is linear with respect to the forces, 
a relation which is not always true. 

A good example of a case in which this linear relation does not hold is 
that of chemical reaction. The thermodynamic forces in this case are the 
differences in chemical potentials, which in turn are linear with respect to 
the logarithm of the concentrations. The rate of a reaction, on the other 
hand, is linear with respect to the concentrations. Hence the rate is linear 
with respect to the exponential of the thermodynamic force, and can be 
written ys, L;;X; only if the concentrations approximate equilibrium 

7 


concentrations (Denbigh, /oc. cit.). 

It is of importance to point out that biological systems are not in quasi- 
equilibrium, but in the steady state. Hence they are not characterized by a 
minimization of 6 by virtue of (6). Actually, in the steady state @ can take 
on any value, as long as it is balanced by a compensatory outflow. 

It might be reasonable, from an energetic standpoint, to postulate that 
living systems do minimize 6. To this author’s knowledge no evidence 
exists to support this attractive hypothesis. 

A fair amount of work has been done with respect to the kinetics of open 
and steady state systems and certain useful things have been learned. 
These have generally been of the nature of finding that certain properties 
of living systems are possible in open systems, rather than that a particular 
property follows necessarily from the fact of an open, steady state system. 
For example, Burton (Joc. cit.) has shown that the “overshoot” phenome- 
non is possible, given certain simple kinetic models of the steady state. 
M. Moore (1949) has shown that autocatalytic systems can give rise to 
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oscillations, etc. Yet it must be pointed out that none of these models lays 
claim to being a replica of an actual living system; hence few conclusions 
can be drawn with respect to the behavior of living systems from the above 
studies. 

There is one further point worth mentioning. L. von Bertalanffy (loc. 
cit.) has defined a concept, “equifinality,” and ascribed it as a general 
property to steady state systems: 

Vital phenomena show a different behavior. Here to a wide extent, the final 
state may be reached from different initial conditions and in different ways. . . . 
In an open reaction system, irrespective of the concentrations in the beginning 
or at any other time, the steady state values will always be the same, being 
determined only by the constants of reaction and of the inflow and outflow. 


As the above statement stands, it cannot be completely accepted. There 
is no doubt that many kinetic models of steady state systems will show 
the property of equifinality. But there exist also steady state systems of 
more complicated structure which will exhibit multiple steady states, with 
various degrees of stability. For example, Denbigh, Hicks, and Page 
(1948) have shown that a fairly simple autocatalytic system can give rise 
to two steady states, one stable and the other unstable. It can also be 
shown very easily that if we assume a system of first-order reactions in 
which all intermediates diffuse, nonlinearities arise, which, under some 
conditions, give rise to multiple steady states. This deserves to be stressed 
for two reasons: 

1) Certain phenomena of seemingly non-genetic changes could be ex- 
plained on the basis of such multiple steady states, whereby environmen- 
tal changes, if of sufficient magnitude, could give rise to shifts from one 
steady state to another and be inheritable. 

2) Most work done so far on the kinetics of the steady state has been 
done on fairly simple first-order systems. Yet there is little doubt that the 
kinetics of living systems is more complicated, probably involving spatial 
parameters, phase boundary discontinuities, and mainly second-order 
reactions. We must therefore be very careful not to generalize too much 
from the properties of the simple models. 

To conclude, then: The kinetic and thermodynamic specifications of 
open and steady state systems have been mainly of use in defining the 
boundaries of what is possible, rather than telling us what actually hap- 
pens. It would seem therefore that more significant information can only 
be gained by a detailed analysis of concrete models. 


_The author wishes ‘to express his gratitude to Dr. John Z. Hearon for 
his many helpful discussions. 
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Rashevsky’s treatment of the flow of an incompressible viscous fluid in an elastic distensible 
tube is applied to the same problem, except that the wall of the tube is assumed to be elastomer- 
ic. As a result the velocity of propagation is obtained in terms of the elastomeric constants of 
the wall, the thickness and density of the wall, the viscosity of the fluid, and the radius of the 
tube. 


The problem of the flow of an incompressible fluid through a distensible 
tube is of great importance in the theory of blood circulation. Several in- 
vestigators have studied the propagation of pressure waves in such a sys- 
tem. In most cases in the literature (Young, 1808; Moens, 1878; Korte- 
weg, 1878; Witzig, 1914; Karreman, 1952) Hooke’s law was assumed to 
hold for the wall of the tube. However, as pointed out by A. L. King (1946, 
1950) this assumption does not hold for the walls of large arteries. This 
author, for the same reason, criticizes the work of N. Rashevsky (1945a, 
b), which showed that the vibrations of the wall of a distensible tube, 
through which a viscous, incompressible fluid flows, may have fre- 
quencies in the audible range. J. J. Hermans (1952) treated the problem 
of the flow of an incompressible viscous fluid through a distensible tube, 
but he neglected the inertia of the fluid. It is the purpose of the present 
paper to show that the method of Rashevsky (1945a) can easily be ex- 
tended to the general case of a viscous fluid flowing through a distensible 
tube with an elastomeric wall, taking into account the inertia of the wall 
as well as the viscosity of the fluid. The influence of the last two factors 
will be shown, and it will be seen that they can be neglected for large 
arteries whose walls are not too thick. In this case a formula which is es- 
sentially the same as that of King (1947) will be obtained. 
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According to King (1946), the relation between the internal pressure p 
and the radius r of a distensible tube with an elastomeric wall is given by 


= oh hn] (29) eG) @ ap c 


In this equation fp is the external pressure, /) and 7p the thickness and ra- 
dius of the unstretched wall, and 6 an age-dependent factor which repre- 
sents the amount of coiling and twisting of the peripheral molecular 
chains. The inverse of the Langevin function, 


L(x) =coth e~s, (2) 


ses. 

We will now apply the dependence of # on 7, as given by (1), to Rashev- 
sky’s work (1945a) to replace his equation (3). Then we have instead of 
his equations (5) the following 


Open 
and 


B | nyt ee) ez 2)", fe 


2179 3 


in which p;', p, are the pressures and r,", r; the radii of two consecutive 
segments of the tube in the stationary state (see Rashevsky’s Fig. 1). As 
in Rashevsky’s treatment, it is supposed that r* varies so little along the 
axis that at each point (1) can be applied. 

The instantaneous pressures at those two points are again denoted by 
pi, px, the excess pressures by (Ap);, (Ap)2 so that 


pi= pi + (AP). (5) 
po= pr + (AD) >. (6) 
Similar to Rashevsky’s equations (7) and (8) we have 
n=nt+A (7) 
ro= 12 +A. (8) 
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The equations of motion are now 


eyeps 
Sp renf nn Bel oy)" 0 


2a rile pi 


and 


on, G 
rot nerd n-ne)" (29 ]ao 


with the same notations as in Rashevsky (1945a), except for the use of 
e instead of 5. Simplifying equations (9) and (10), substituting the results 
obtained from (3), (4), (5), and (6) by eliminating p; and p, we have, 
also using (7) and (8), 


€ pi = (Ap)i+ me ~ Ry ae * C5 ws) 


ae 


ra(s2) a 
int B 5 ne 
22 cys 
and 

: nt (e) 
ep: oo = pest (apr eae) oan (12) 


ie oan fe (6 fs) =p: 2)": 

r i 

If we neglect the inertia of the wall, (11) and (12) reduce to equation (10) 
of King (1947), who treated this case in a different manner: 


i (5 *) a2 ee i> (87 “) 


_ lobo ipl 
ap att ( eat a Lat (ane Etat Ges (13) 
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Using (7) and (8) 
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may be expanded about r;. Keeping only the first two terms, we have 


(King, 1947) 
(eta (o\r40%), ona a9 


with 
ux sinh? uz 
= 9 ob Se = 5 
= B nha * — a? Delis 
and 
Peay pad (a ar), Pa be ei 
To 


Substitution of (7), (8), and (14) into (11) and (12), keeping only terms up 
to A,/r; for A; < r; (4 = 1, 2), leads to 


eo Sat = eae ei" ma te -1) 
(17) 
+3 LM Ai 
ry 
and 
Se aCe 
(18) 


3/2 
a3 7) A 
2 


If we consider with Rashevsky (1945a) only deformations which preserve 
the constancy of the volume, we have 

Ar= —As.« (19) 
If we subtract (18) from ate and use (19), assuming that r;" and 7; are 
close to their average value r*, we find 


ses sy (Le aaa 1(8 2") 


2ep, At = (Ap) 1— (Ap) 2— dr hfe 


(20) 


x (2B = -1)+3(4 1 
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This equation is the equivalent of Rashevsky’s equation (12). Elimination 
of (Ap): — (Ap)s from (20) and Rashevsky’s next equation (1945a, p. 30), 


_ pP day | 8nl? dA 
rr or rh = (Ap) 1— (Ap) 2+ 7S ore 


(21) 


(with r* instead of his R*), which is an approximate form of the Navier- 
Stokes equation, 7 being the coefficient of viscosity, leads to 


(rent) fii Sah dar, tte ((raynt (7) 
dt? r*8 dit Wi iyere LE (8) 
(22) 
r* ro \3/2 
This equation allows a solution of the form 
Ay = el (nl?) /r#*(2¢p,7*+01*)1¢ (4 cos w t+B sin wt) , (23) 
in which 
ans loPo (2 inl (8 Z) 
r ancacttncar| = E28) (25 a) 
(24) 


1/2 
Re eye (a in ae 
aN Toa 2epir*8-+ pl?r *2 , 


Since the wave length is 2/, we have for the velocity of propagation of 
the disturbance 
J 2 : (5) 
T 
The velocity of propagation is given by (24) and (25) in terms of the 
elastomeric constants, the thickness and density of the wall, the radius of 
the tube, and the viscosity of the fluid. 

From (24) we see that for the large arteries the influence of the viscosity 
on w and, therefore, on 2 is very slight, since for those vessels 7» and r* are 
of the order of 1 cm., é9 and e of the order of 10-1 cm., whereas fo is of the 
order of 10° dyne cm.~, p and p1 of the order of 1 gm. cm.~%, / of the order 
of 1 cm., and 8 of the order of 1. Therefore the second term within the 
outer bracket in (24) is very small compared with the first one, which 
shows that the influence of 7 on w and, therefore, on v can be neglected. 
The damping due to the viscosity (n = 0.02 gem.* sec.) as given by 
(23) may play a role for slightly smaller arteries, for which r* is of the 
order of 10-1 cm. The influence of the inertia of the wall depends on the 
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order of the term 2epr* compared with p/? and therefore approximately 


on that of 2er* compared with /’. mae 
For a non-viscous incompressible fluid 7 = 0, and if the inertia of the 


wall can be neglected (24) and (25) reduce to 


Oe OD) (apts) a8 a \y (26) 


The dilatation r*/ro is then given by (3) or (4) with ph, = p = 'p* and 


r=r=r* 
*_p = oho s(n a ~) -(4 s) f (27) 
‘ TL ® 
If we put r*/ro = x, x is found by solution of 


Ore ae), ae 


which yields x as a function of p*. 

For example, in the case of a 50-year-old person we have (King, 1946, 
1950): 6 = 0.52, ro = 6.6 mm., e9 = 1.12 mm., and graphical solution 
of (28) for x leads to Table I. 


TABLE I 
(o*—p0) in mm. Hg x (/*—}o) in mm. Hg x 
10 1.07 70 1.42 
20 1.14 80 1.46 
30 he20 90 1.49 
40 12% 100 1.52 
50 ESS 110 255 
60 eos 120 1.58 


From Table I we see that for an intra-arterial pressure of 70 mm. above 
atmospheric pressure the arterial wall is stretched 42%. 


Elimination of 
Cee 
y 8) 
from (26) and (27) yields 
Ms 1 lopo 2ro p* = Po 3/2 jal 50 
Join teed ED I Po par ype te r* {( 25 145 & ; sabe 


This formula is the same as that of King (1947), except for the numerical 
factor in front of the first brace. 
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If ro/ly is constant, which might be a reasonable approximation as a 
first step, x according to (28) is independent of 79. Then we obtain from 
(29) 

A 
a ANTE! (30) 


in which A depends on é, po, p*, p, and 8. 

This investigation was supported by a research grant H-627(C2) from 
The National Heart Institute, of the National Institutes of Health, Public 
Health Service. 
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A derivation is given of the reflection coefficient of pressure waves in a vessel whose end 
branches into many smaller vessels. This coefficient depends on the number of these smaller 
vessels and their sizes relative to the size of the main vessel. Estimations are made of the order 
of magnitude of the coefficient. Assuming the main vessel to be of the order of size of an artery, 
it is shown that the reflection coefficient has a value close to one for reflections at branchings 
into vessels of arteriolar size. It is pointed out that the result may support the idea that the 
standing waves in the arterial system are due to reflections at the site of the arterioles. 


W. F. Hamilton and P. Dow (1939) discovered the standing pressure 
wave system in the aorta. This discovery implies that the pressure waves 
are reflected in the arterial system (Karreman, 1952). It seems obvious 
that this phenomenon would occur in a system with many discontinuities, 
such as the circulatory system. The importance of this phenomenon was 
also emphasized by A. Apéria (1940) and J. G. Porjé (1946). As was shown 
earlier (Karreman, 1954), by using certain assumptions and approxima- 
tions, * the velocity v of propagation of pressure waves in large tubes may 
be given by 


A 
eye (1) 
where r represents the radius of the tube. 
Consider the system sketched in Figure 1. Continuity of the pressure at 
the site of branching is expressed by 
Psa=?Pc (2) 
and continuity of flow by 
TR20p =NTL* VC , (3) 
where v, and v¢ represent the velocity of the fluid at B and C. If  repre- 
sents the pressure of a pressure wave traveling from A to B and its reflec- 


* It was noted in this paper that this result is valid if the influence of the viscosity of the 
fluid and the inertia of the wall are negligible. Crude estimations show that this is the case 
down to vessels with radius of the order of 20 yp. 
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tion at the site of branching, then we have for the pressure and velocity of 
the fluid in AB (Karreman, 1952) 


por(i-ri(49 : 


-2r(-3)-Ls(49). 6 


If we suppose that there is no reflected wave in CD, we have (Karreman, 
1952) for the pressure and velocity of the fluid in CD 


p=e(t-4) (6) 


1 x = 
where ¢ and c’ are the velocities of propagation of the pressure waves in 
AB and CD respectively. 


and 


2R ii ae ss n 
La Bie YO vee 
| * 2A 
[i ane Uehe [eaaEE Y 
A B 
Ficure 1 


If we put the origin at the site of branching, we obtain from (4), (5), (6), 
and (7), by substituting 0 for x, 


ba=F(t) +7 (0) (8) 
1 

tet) — FI (9) 

c= g(t) (10) 
1 

oe Aaa San (11) 


Substitution of (8), (9), (10), and (11) into (2) and (3) leads to 
F(t) +f (t) = g(t) (12) 
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and 
R?2 
(FO —fOj=n5 ew. (13) 


If, as before (Karreman, 1952), we define the reflection coefficient n by 
f (2) 


we obtain from (12) and (13) 
{284) fake 
R? cf 
= =F (15) 
1 +n 5s ; 
Application of equation (1) leads to 
necte 
c VR A/fr\2 
sori (16) 
Vr 
Substitution of (16) into (15) yields 
= Ai \R 17 
if 14 A Y 5/2 ( ) 
"aa z) 


If we assume that Hooke’s law is applicable to the wall and use the 

approximate formula (Karreman, 1952) for the velocity of propagation 
eis 

fava pr’ 


in which £ is the elasticity modulus, 6 the thickness of the wall, and p the 
density of the fluid, we have 


Similarly 


Substituting these expressions for A and A, in (17), we obtain 


1 ey) 1G)" is (18) 
Fy6 aa (a) 
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From this formula we see that if mr/?>> R°? and E = &,, 6 = &, or, if 
Ey, <E and nr’/? = R°/2, § = &, we have 7 = —1, corresponding to a 
phase change of 180° at the (complete) reflection. This is the case if there 
is a localized increase in the total cross-sectional area or a localized in- 
crease in elasticity of the wall. For nr? << R°? (E= &, 6 = &) or 
E, > E (nr°/? = R°/2, § = 6), corresponding to a localized decrease in the 
total cross-sectional area or a localized decrease in the elasticity, we have 
n = +1, corresponding to a complete reflection without a phase change. 
For reasons of simplicity, let us suppose that 


A = A if". (19) 
Substitution of (19) into (17) gives us 


es fe 
i+n(g) 


H. D. Green (1950) has compiled a table of the sizes and numbers of 
vessels of a certain kind in the dog. We will take for the average size of the 
vessel in front of the site of branching in Figure 1 that of an artery witha 
radius of 0.3 cm. Therefore we assume R = 0.3 cm. If we further assume 
that the branches are of the size of the “secondary branches,” we find 
from Green’s table that r = 0.03 cm. and m = 1800. Substituting these 
values in (20), we find 


U] (20) 


pee 0270 « (21) 


Assuming that the branches are “tertiary branches,’’ we have r = 0.007 
cm. and x = 76,000; therefore, for this case 


ee 0. Es (22) 
For terminal arteries r = 0.0025 cm. and m = 105, so that 
n= —0.73, (23) 
In the case of ‘terminal branches” r = 0.0015 and = 1.3 X 107, so that 
n= —0.92., (24) 
Finally, for arterioles, r = 0.001 cm. and = 4 X 10", which leads to 
n=—0.92. (25) 


Porjé (Joc. cit.) has shown that more than 98% of the energy of the 
pulse is accounted for by the first three components of its Fourier analysis. 
In most cases, the fundamental component accounts for about 80% of the 
total energy of the pulse. For this reason we simply study the pressure pat- 
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tern due to one component wave and its reflection. Suppose, therefore, 
that the incident wave is given by 


p=A cos w (1 *). (26) 
The reflected wave will then be given by 
p=nA cos w (14+). (27) 


The total pressure due to superposition of the two waves (26) and (27) will 
then be determined by 


P= A}cos w(1—=)+n cos (/+*)} | (28) 


= A (14) cos wf cose + A (1— 7) sin of sin w =. 


As was pointed out formerly by this author (1952), the pressure pattern 
(28) never corresponds to a standing wave, unless one of the two terms of 
the right-hand side of (28) vanishes, which happens only if 7 = +1 or 
— 1, corresponding to complete reflection just in or just out of phase (i.e., 
with a phase difference of 0° or 180°). Now we see from (21) through (25) 
that 7 is closest to —1 for the case of vessels approximately the size of 
arterioles. In this case, the amplitude of the second term in (28) is 1.92A, 
whereas that of the first term is 0.08A, so that this second amplitude is 
only a little more than 4% of the first one and may not be detectable. 
This analysis might then give some support to the claim of Hamilton and 
Dow (loc. cit.) that the standing wave pattern is due to reflections at the 
subterminal branches. 

Of course, we realize that the model treated here is a physical more than 
a biological one. In the latter, the branches between an artery and the 
arterioles should also be taken into account. In addition, even in the 
physical model, the pattern due to the reflection at the smaller vessels 
should be investigated in much more detail. 

Finally we would like to point out that the damping due to the viscosity 
can be approximately represented by a multiplicative factor e740) /(oR)t 
(Karreman, 1954). Taking this into account, we have, instead of (28), 


pA e~ (nfo?) e} (iy) cos'ot cos @ a+ (A =a) sinotsin ={. (29) 


This investigation was supported by a research grant H-627(C2) from 
The National Heart Institute, of the National Institutes of Health, Public 


Health Service. 
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The paper is a second step toward a biomathematical theory of the rates of spread of new 
nonconformist ideas or behaviors in a society. It is intended as a preliminary and purely theo- 
retical study of a very oversimplified case. An equation which determines the distribution 
function of the tendencies toward conformist and nonconformist behaviors is set up under a 
number of oversimplified assumptions, and a solution by successive approximations is indi- 
cated. The expression for the first approximation is given, and an estimate of the order of mag- 
nitude of the rates of changes is made. In conclusion an outline is given for further improvement 
of the theory. 


In a previous paper (Rashevsky, 1953) we discussed a biosocial mecha- 
nism which may account for the gradual spread of new ideas or forms of 
behavior through human society. In particular we discussed the gradual 
change of society from “‘arational” behavior (Rashevsky, 1951), which is 
based principally on an acceptance of a set of beliefs or prejudices, to a 
“ational” behavior, based on a critical and logical examination of any 
action or opinion. The gross features of development of humanity from 
the stage at which it was some 8,000 to 10,000 years B.c. to our days do 
show just such a slow shift from arational to rational behavior. 

The mechanism discussed is based on a general assumption, which 
hardly will be doubted by any biologist; namely, we consider that the be- 
havior of every individual is determined both by his innate tendencies, 
and by the educational effects of his social environment. The natural 
spread of the innate tendencies results in itself in a non-uniformity of be- 
havior of a social group. A certain fraction of individuals will “naturally” 
behave appreciably differently from the average majority even in the ab- 
sence of any social environment. This fraction of “nonconformists,” how- 
ever, affects the social milieu in which new-born individuals are raised, 
and therefore an individual of a new generation, who has the same innate 
tendencies as another individual of the old generation, will in general be- 
have differently from that other individual. This results in a change of the 
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fraction of the nonconformists, and this in its turn results in a further 
change of the behavior distribution of the next generation. 

We have derived (1953, p. 222) an approximate expression for the 
rate of change of the fraction of nonconformists for the case of two mutu- 
ally exclusive behaviors (such an arational and rational). Implicit in that 
derivation is the assumption that the distribution function for the total 
tendency toward one or other of the behaviors, that is, the tendency 
determined by both innate and environmental factors, merely shifts 
toward one of the behaviors, without changing its shape. This assump- 
tion is most certainly unrealistic, and should be discarded in the next step 
of our investigation. Setting up the equations which govern the time 
course of the distribution function under most general assumptions does 
not present too great difficulties. The nonlinear functional equations thus 
attained are, however, of a kind which hardly lends itself to treatment by 
known mathematical methods. The study of some rather oversimplified 
and unrealistic cases, which are amenable to analytical treatment, is 
therefore indicated at this stage. The study of such purely theoretical 
cases will give a general orientation in this field and may lead to expres- 
sions which may be used as zero order approximations to the more realistic 
but also much more difficult situations. All that could be expected from 
such a preliminary theoretical study by way of practical application are 
correct orders of magnitude when plausible values of the parameters 
are assumed. 

Let ¢ denote the difference between the central excitations « and e, 
which determine the drives toward the two mutually exclusive behaviors 
R, and Ry. When an individual has the value ¢ = 0, he is indifferent in the 
choice of R; and R. A large positive ¢ means a strong preference to Ri, a 
large negative ¢ a strong preference to Rp (Rashevsky, 1951, Chaps. ii, 
K1151953): 

Let the innate distribution function of ¢ among newborn individuals 
not yet influenced by the social environment be 


No(¢) do (1) 
with 
+0 
SJ. Nolo) dg =1. (2) 
In the total population of individuals of different ages the value of ¢ of 
each individual is in general determined both by innate factors and by the 
effect of education. Let 


N (4, ¢) (3) 
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denote the distribution function of ¢ for the whole population at the time 
t. For all values of t we have 


+c 
ie N (¢, t) dg@=1. (4) 


We now shall make the following highly oversimplified and unrealistic 
assumption. 

Let m denote the average number of persons with whom a newborn in- 
dividual comes in contact and who influence his #. We shall consider that 
this influence lasts a period of time during which the newborn individual 
grows up, and that after that the ¢ of the individuals does not change 
any more. Further we shall consider the case in which N (4, #) changes so 
little during the period of “growing up,” that on the time scale chosen this 
period may be considered as infinitesimal. 

Furthermore we shall make the following assumptions about the effect 
of influence of an adult with a value ¢” upon a newborn with value ¢’: 


a) If ¢’ > ¢” there is no effect. 
b) If ¢’ < ¢” then the value ¢’ changes to 


Pai? — oe), (5) 
where a is a constant, and 
Demi” (6) 


The value of ¢” of the adult remains unchanged. 

Assumption a) appears to be so artificial as to be almost nonsensical. 
This is not so, however, when we consider the type of behaviors R; and R2 
in which we are interested. R; stands for rational behavior, R2 for aration- 
al. A young individual with a large positive ¢’ is naturally skeptical of 
the accepted type of behavior which he is being taught. If he comes into 
contact with an adult whose ¢” < ¢’, that is, one who accepts without 
criticism the behavior accepted by society, the young individual is not 
likely to be in sympathy with him, and is not likely to be influenced by - 
him. On the other hand, if the young individual meets an adult with 
¢” > ¢’, that is, one, who, like himself, is a skeptic, his skepticism is likely 
to be enhanced. The choice of (5) for the increase of ¢’ in this case is arbi- 
trary, and merely represents one of the simplest possible assumptions. A 
different, more plausible assumption will be further discussed below. 

The biosocial picture outlined above should actually lead, for the case 
¢’ > ¢”, to a change from ¢’ 


¢’ — a’ (¢’— $”) 
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ith 
> Osa <a ie (7) 
Assumption a) represents a limiting case, when a’ becomes negligible com- 


pared to a. 
Instead of speaking of “individuals with a value ¢” we shall speak of 


‘Individuals ¢.” , 
The probability that an individual ¢’ has one adult ¢” among those 
with whom he comes into contact is given by 


mN ($"’, t) dg”. (8) 

The probability of having one or more adults ¢’’ among those he meets is 

1— [1—N (¢”, t) do’’"] ™ = mN (¢”, t) do’’ + infinitesimals of higher order , 
and need therefore not be considered. 

Let the relative birth rate u and death rate v for the society be constant 


and independent of ¢. This is the simplest possible assumption. If 3(¢) 
denotes the total population at the time #, and if 


ae (9) 
then 
oe = (u—») R(t) = «N(); soe 
so that 
MN (2) =MNoex*. oe 


The number of individuals ¢’ born during the interval #, ¢+ dt is 


equal to 
uN (t) No ($’) do’dt. (12) 


During this time dé, which, according to what we said on page 119, is 
also their “growing up” time, they will meet, according to (8), 
um No (o’) N (o”, t) db’, do’ dt (13) 
adults $””. 
Provided $” > ¢, the value ¢’ of the newborn individuals, whose num- 
ber is given by (12), will shift to 


o= ob ale 6). (14) 
Putting 
B=1—-a<1, (15) 


we have from (14) 


ge” =——* ; dg” =— dg. (16) 
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Because of (16), expression (13) may be written, omitting, for brevity, 
the argument ¢ in J, 


ME M(t) No(6!) N (e—**) do'dodt. (17) 


The total number of individuals who shift during dé into the interval ¢, 
¢ + d¢ is obtained by integrating (17) with respect to ¢’ from — ~ to 
¢, since, according to assumptions a) and b), p. 119, only upward shifts 
occur. Hence this total number is given by 


C.dodt=™* 7 (1) dgat [° No (6!) N (Ssh) dd’, (18) 
a —o a 


Because of (12), the number of individuals ¢ born during di is given by 


Crdgdt = pM (t) No (o) dod. (19) 
The loss of individuals ¢ during dt through death is 
Ladgdt =vM (1) No($) dodt . (20) 


In addition to the above contributions and losses, there is a loss of new- 
born individuals ¢ due to the contacts with adults with ¢” > @. This loss 
is equal to the total number of encounters between newborn ¢ and 
¢” > ¢, or, because of (13), to 


L,dodt = mp (t) No(s) deat ih “N (6!) do”. (21) 


The total rate of change of all individuals 2(#)N(¢, t)d@ who have the 


value ¢ is 
dt (t) N (¢, 4) 
dt 


or, because of (9), (18), (19), (20), and (21), 
ets? =n f No($!) N (S@—* Be ag" 


=C,+0,-—L—-fL,, (22) 


(23) 
+ M(t) No() — muN (4) No (4) hi N (¢") do”. 


We have 
dN ee t) =N (1) oe t) +N (¢, 2) ae 2 


? 


or, because of (10), 
PEO AS = RAY) [eet vig, |. (24) 
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Put 
ms ph eee (25) 


Introducing (24) and (25) into (23) we find 


SHAG DE {Ne (° =") N (x; 1) dx 


(26) 
bxtiden's () fN(e, EAC AC hONE 
d 


which is the functional equation for the determination of the unknown 
function V(¢, t), when NV,(¢) is given. 

A different equation is obtained if we consider the following more real- 
istic picture, instead of expression (5): 

When a newborn ¢’ meets an adult ¢” > ¢’, there is a probability 


b (A; ¢’, o’’) (27) 


that his ¢’ will be increased by A (> 0), which is positive, but less than 
go” — ¢’. We have 


fO “pda= C (28) 
0 


The requirement that A < ¢” — @¢’ is essentially the mathematical 
expression of the biblical saying that the disciple is not better than the 
teacher. If, as seems to be actually the case, a “‘disciple” may well exceed 
the “teacher,” then the restriction A < ¢” — ¢’ may be dropped, and 
the integration in (28) taken from 0 to +. 

Of all newborn individuals ¢’ meeting adults ¢”, as given by (13), the 
number 


mp M(t) No(o’, t) N (¢"’, t) bp (6 —¢'; 6’, 6”) do’dg’'dodt (29) 
will shift to ¢. Hence, instead of (18), we now have 
¢ pte 
Com mun () [" f No(o') N (6's) (6-9; 8", 6") do's”. (30) 
Instead of (26) we now obtain 
dN (¢, t) ¢ ER / ad , / v7 
a ff Noo") N (6, 1» (6-4; 6°, 6”) dg'de” 


= (31) 
= miko One N (x, t) dx. 


We shall discuss here a method of solving equation (26) by successive 
approximations. The same method can be used for equation (31). 
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Equation (26) may be 


a mu f [5 No (25 **) — wo) ] N(x, t) dx. (32) 


Differentiating 7 times with respect to ¢, we find 


N® (¢, 0) = oN = mu f LE No (25**) — woe) | 


dt* 
m (33) 
ye EIN (a, 


eet ee) O52) 


If No(#) and the function (4, 0) are analytical, then all the derivatives 
N(¢, 0) can be successively calculated, and we may express NV(, #) as 
a Taylor series in #: 


N (4, ft) = N (gd, 0) + NM (g, 0) t+ 9N (o, O)P+.... (34) 


We shall consider the case where Vo(¢) is a normal distribution. Setting 
a zero point for the time in such a problem as this is always rather artifi- 
cial. It is, however, not unnatural to assume that in the earliest days of 
our civilization, when population densities were small and social contacts 
very few, N(¢, ¢) was essentially the same as Vo(¢). Hence we shall put 


N (¢, 0) = No (¢) -+ g $e oe) e-l@te)*/e7], (35) 


We thus consider in general the possibility of an initial bias ¢o z 0. 
Introducing (35) into (32) we find: 


“g (See) (=) dg 
me (#48) f°" (CY) on 


¢=o+ Boo (37) 


we may write the integrand of the first integral in (36) thus: 
o— ax «+ go 
ca ghaee). We) | 38) 


B 
1 art B se 


N’ (¢, 0) Say aie 
(36) 


Putting 


Put 
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and 
_ B(1—2a)do— ae 
i= a2-+ 6? 
Vat Bo=b. 


Now (38) may be written 


1 (*) 
BN OTy)/ PJ ——— ]}. 
o2 ar : 8 oat 


Hence (36) may be written 


=e (29) f(a 


Making the substitution 
er Oi _ 


O71 y 


in the first integral of (43), and the substitution 


e+ Go _ . 
o 


in the second, we find 


, _ Mio, o+ do tz 
W'($, 0) =H g (Se) fal) dy 


me co) +e) Pow 
at. (SE ds. 
3 5 ( o as ( 2) : 


Introducing the notation 


G(u) =f g(y) dy; [G(@) =1;G (0) = 4) 
we find 
+o 
fs) dy=1-Gu). 

The relation between G(u) and the function 

2 a ee 

®(u) = BF i di, dx, 

which is available in tables, is 


G(u) =3[14+(u)]. 


(40) 


(41) 


(42) 


(43) 


(44) 


(45) 


(46) 


(47) 


(48) 


(49) 


(50) 
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Hence, from (48) and (50) 


J sQ)dy=4 1-0 (u)]. (51) 
Introducing (51) into (46) ,and remembering (39) and (40), we find 
! = mi d+ go o+¢ 
WG, 0) = 7 aes & (2) {1a (1-20) 25%] 


Be , (£8) [1-2 (224), 


Because of (15) 8 can be expressed in terms of a. Using (33) we now can 
compute WV’’(¢, 0), and so on. For arbitrary values of a between zero and 
one, expressions for the higher derivatives cannot be obtained in closed 
form. 

What interests us here, however, is not NV(@, #) itself, but the fraction 
of nonconformists, which is determined by N(¢, #). That fraction y is 
given (Rashevsky, 1951, Chap. xii) by 


+o 
y= f_ N(@,) 1-Pi) 1 dg, (53) 


where P;(¢) is determined by the distribution function p(£)dé of the 
physiological fluctuations within an individual, which result in the vari- 
ability of his behavior (Landahl, 1938; Rashevsky, 1948, Chap. xxxiv; 
1951, Chaps. ii, xii). A closed expression for y has been given by H. G. 
Landau (1950) for the case in which both V(4, #) and p(é) are normal 
distributions. Since in the present instance NV(¢, ¢) does not remain nor- 
mal, even if it is so at ¢ = 0, therefore, we may again limit ourselves to a 
special case, namely, where the standard deviation of p(£) is very small, 
in other words, where the inner fluctuations are negligible. In that case an 
individual with ¢ > 0 always exhibits behavior Ri, one with ¢ < 0 always 
exhibits Ro. Therefore, if R; stands for nonconformist behavior, the frac- 
tion y is simply equal to 


= [ N(@,t) dé. 54 
y= fo NG, 0 de (54) 
Consider in particular the initial distribution (35). Then the initial value 
yo of y is given by 
wl? of B28 
y= fg (22) ae, (55) 


or, using (51), 


vo=4[1-2(2)]. (56) 
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We have previously (Rashevsky, 1953) estimated yo at 10-%. Using 
tables for ® (Peirce, 1929, p. 116), we find 
BES I (57) 
oO 


Expression (52) represents the increment of V (¢, ¢) at ¢ = O per unit 
time. Therefore 


GY | [ Nea ne (58) 
Tl, = J, N’ (0 de 
represents the increment of y during that interval. 

In general we have 


diy ~d' N (9; #) 
ae =e ee LEM Ata ears ©?) 


We therefore can compute successively all the derivatives d‘y/dt* for 
i = 0 and use a Taylor expansion for y, just as in (34). 

For a = 0 we have N’(¢, 0) = 0, as should be the case. For a = 1, 
which gives the fastest possible rate of change of V(¢, #), we find, because 
of (41), 


N’(g, 0) = FEY (PFS) 1 +9 (SE) 


Grab eae it 


The first expression in brackets is very close to 2 for ¢ > 0. 

For a = 1 all the derivatives N° (¢, 0) can be expressed in closed form, 
for their computation reduces to evaluation of integrals of g(x) and of 
products g(x)@(«). The successive derivatives do not show, however, any 
regularity in their form, which would permit a guess at a general expres- 
sion. An elaboration of this case is hardly worth while, because the as- 
sumption a = 1 is very implausible. Inasmuch as nothing is known about 
a, we cannot estimate the order of magnitude of the time ¢* during which y 
increases from ‘yo to the order of magnitude of unity. Very crude numerical 
estimates indicate that with a ~ 0.1, which is as plausible a value as any, 
the value é* is of the order of 10 years. This of course does not mean any- 
thing, except that correct values of #* may be obtained with not implau- 
sible values of a. A further elaboration of the theory is necessary before 
more definite conclusions may be made. 

As we said, it is highly unlikely that a would be of the order of one. Such 
an order of magnitude would imply that almost every child who comes in 
contact with one nonconformist adult would grow up into a nonconform- 
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ist. If out of m adults who influence the child there is one nonconformist, 
a is more likely to be of the order of 1/m. In some cases a nonconformist is 
likely to be more influential than a conformist, and contribute more than 
1/mth of the total influence. If his effectiveness is that of k < m individu- 
als, then a may be approximately given by 


Se h 
 mEh—1? 


a (61) 
which reduces to 1/m for h = 1. 

The greatest oversimplification is the assumption that N/(¢, £) changes 
very little during the “growing up” or “formative” period of an individu- 
al. Actually an adult individual may also change his ¢ under the influence 
of others, though the probability of such a change decreases with age. 

As a next step in the development of this theory, which is bound to lead 
to mathematical problems of great difficulty and possible mathematical 
interest, we envisage the following: 

A distribution function V(q, 7, #) should be considered, which gives the 
fraction of individuals of given ¢ and of age 7 at the time ¢. Contact of an 
individual ¢’, 7’ with an individual @”, 7’ results in a finite rate of change 
of ¢’, so that during the interval di, the value of ¢’ changes to 


1, ag’ 
g ss gs dt, (62) 


where d@’/dt is a function of ¢’, ¢”, 7’, 7’, and m. Thus 


i hoe Tm) (63) 
dt 
The relative birth rate v should be considered as a function of the distri- 


bution of ages 
+a 
M(r,) =f N(o,7,0 46, (64) 


and the death rate p should be considered a function (7) of the age, and 
possibly even a function y(r, ¢) of the age and of ¢. The dependence on 
may be quite complex. On one hand, a large positive ¢, which charac- 
terizes an individual with strongly rational behavior, may help the sur- 
vival; on the other hand, in a preponderantly arational society, rational 
nonconformists with a large positive ¢ might be persecuted and even ex- 
terminated, as has actually happened and is still happening. 

If, with all these considerations in mind, we establish the equation 
which governs the variation of N(¢, 7, t) with time, we shall have come 
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closer to real situations, but also closer to almost insuperable mathemati- 


cal difficulties. 
Problems of history may still turn out to be an inspiration to mathe- 


maticians, as problems of physics have been, and as problems of biology 
are bound to become. 


The author is indebted to Professor George Karreman for a critical dis- 
cussion of this paper. 
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